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A NUMERICAL  SOLUTION  OF 
SUPERSONIC  AND  HYPERSONIC  VISCOUS 
FLOW  FIELDS  AROUND  THIN 
PLANAR  DELTA  WINOS 
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A numerical  technique  was  used  to  compute  the  supersonic  and 
hypersonic,  viscous  flow  fields  around  thin  planar  delta  wings. 

These  solutions  were  obtained  by  solving  the  Navier-Stokes  equations 
subject  to  a conical  approximation.  The  integration  technique  used 
was  the  second-order  accurate,  finite-difference  scheme  by  MacCormack. 
This  numerical  integration  was  performed  on  a constant  step  size  array 
generated  by  a conical  coordinate  transformation.  Solutions  were 
obtained  for  the  upper-only,  lower-only,  and  total  flow  fields  around 
delta  wings  with  supersonic  leading  edges.  These  solutions  span  a 
Mach  number  range  of  2.94  to  10.17,  a local  Reynolds  number  range  of 
3.345  x lo"'  to  5.0  x 10^,  and  various  angles  of  attack  from  -151  to 
+15°.  Numerical  oscillations,  due  to  shock  capturing,  were  reduced 
by  applying  normal  stress  damping  and  a fourth-order  density  damping 
terra  to  the  finite-difference  equations.  A stability  criteria  was 
developed  and  used  which  accounted  for  both  the  viscous  and  tnviscid 
flow  regions.  Good  agreement  was  obtained  between  the  numerical 
results  and  the  experimental  flow  field  data  by  Cross,  Spur 1 in,  and 
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Banntnk.  The  shock-induced  vortex  within  the  viscous  region  and  the 
hypersonic  viscous  "bubble"  on  top  of  the  boundary  layer  were  computed, 
for  the  first  time.  A unique  examination  was  made  of  the  vortical 
singularities  in  the  conical  cross-flow  plane  of  the  delta  wing.  This 
investigation  demonstrated  the  feasibility  of  applying  the  conical 
approximation  to  the  Navier-Stokes  equations  In  order  to  solve  flow 
fields  around  thin  delta  wings. 
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A NUMERICAL  SOLUTION  OF 
SUPERSONIC  ANO  HYPERSONIC  VISCOUS 
FLOW  FIELDS  AROUND  THIN 
PLANAR  DELTA  WINGS 

I.  Introduction 

Aerodynamicists  have  been  investigating  the  supersonic  and  hyper- 
sonic flow  fields  around  delta  wings  for  a number  of  years.  These  flow 
fields  are  essentially  conical  at  high  Reynolds  numbers  and  thus  are  of 
great  interest  both  theoretically  and  experimentally.  Initial  efforts 
at  studying  this  problem  were  directed  towards  analyzing  the  lower  sur- 
face or  compression  side  flow  field.  However,  from  recent  experiments, 
it  was  seen  that  the  most  striking  features  of  the  flow  exist  in  the 
upper  surface  or  leeside  flow  field.  At  angle  of  attack,  an  embedded 
shock  is  formed  on  the  expansion  side  of  the  delta  wing.  The  inter- 
action of  this  embedded  shock  wave  with  the  attached  boundary  layer  at 
hypersonic  speeds  generates  two  important  features  of  the  leeside  flow 
field.  These  include  embedded  streamwise  vortices  in  the  boundary  layer 
at  low  angles  of  attack  and  flow  separation  at  large  angles  of  attack, 
both  of  these  features  increase  the  temperature  and  pressure  gradients 
on  the  upper  surface  of  the  wing.  The  purpose  of  this  investigation  is 
to  develop  numerical  solutions  for  this  flow  field  which  are  valid  for 
both  supersonic  and  hypersonic  speeds.  The  delta  wing  which  is  being 
modeled  in  this  study  is  a thin  planar  wing  with  supersonic  loading 
edges.  Solutions  to  this  problem  for  both  the  upper  and  lower  surface 
flow  fields  are  of  great  practical  importance,  since  lifting  re-entry 
vehicles  such  as  the  Space  Shuttle  have  delta  wing  configurations. 
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Description  of  Flow 


For  a comprehensive  understanding  of  this  complex  flow.  It  is  neces- 
sary to  examine  the  flow  fields  on  both  the  compression  and  expansion 
sides  of  a planar  delta  wing  at  subsonic,  supersonic,  and  hypersonic 
speeds.  Particular  attention  is  placed  on  the  leeward  side  of  the  wing, 
since  this  flow  has  a significant  effect  on  the  lift,  drag,  stability, 
and  control  effectiveness  of  the  vehicle. 

Subsonic  Flow  Field.  Several  investigators  (Refs  1-4)  have  examined 
the  flow  characteristics  around  a slender  delta  wing.  At  subsonic  speeds, 
the  flow  around  a thin  delta  wing  with  sharp  leading  edges  changes 
significantly  with  angle  of  attack.  At  small  angles  of  attack,  the  flow 
remains  attached,  but  as  the  angle  of  attack  is  increased  beyond  some 
critical  value,  the  flow  field  becomes  dominated  by  a large  region  of 
separated  flow.  The  principal  features  of  this  separated  flow  are 
illustrated  in  Figure  1.  The  primary  flow  separation  occurs  at  the 
leading  edge  where  a vortex  sheet  is  formed.  This  vortex  sheet  curves 
upward  and  inward  and  eventually  rolls  up  into  conical  spirals  that 
form  concentrated  vortex  cores  above  the  wing  surface.  The  vortex 
causes  a reattachment  of  the  flow  inboard  of  the  leading  edges,  creating 
two  flow  patterns,  a primary  vortex  flow  and  a secondary  vortex  flow. 

The  primary  vortex  flow  consists  of  concentric  spirals  about  the  core, 
the  core  flow,  and  the  flow  over  the  wing  surface.  The  core  flow  is  a 
highly  rotational  flow  field,  in  which  the  velocity  and  pressure  fields 
arc  approximately  axially  symmetric  (Refs  5-11).  The  viscous  effects  and 
compressibility  effects  are  confined  to  a region  near  the  vortex  axis. 

The  spanwise  flow  next  to  the  wing  surface  is  initially  induced  outward 
toward  the  leading  edge.  Shortly  after  the  fluid  passes  underneath  the 
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vortex  core,  the  adverse  pressure  gradient  causes  the  boundary  layer  to 
separate  and  form  a secondary  vortex  Just  Inboard  and  above  the  leading 
edge.  The  circulation  of  this  very  weak  secondary  vortex  Is  opposite  In 
sense  to  the  primary  vortex  and  has  a circulation  of  approximately  live 
percent  of  the  primary  vortex  (Ref  6). 

Supersonic  Flow  Field.  As  the  Mach  number  Ls  Increased,  the  flow 

on  the  leeward  side  of  the  delta  wing  transitions  from  a separated  flow 

to  one  characterized  by  a Prandtl-Meyer  expansion  in  the  cross-flow 

plane.  The  conditions  under  which  this  transition  occurs  depends  on  a 

number  of  parameters,  such  as  leading  edge  shape,  tiirfoil  profile,  and 

angle  of  attack  (Ref  12).  Intuitively,  the  range  of  conditions  over 

which  transition  occurs  for  various  wing  configurations  Is  approximately 

centered  around  M =1.0  (where  M ls  the  Mach  number  normal  to  the 
n , n 

leading  edge).  For  thin  wings,  the  transition  occurs  well  before  the 
leading  edge  becomes  supersonic.  On  thick  wings,  at  relatively  low 
angles  of  attack,  a strong  detached  shock  wave  forms  in  front  of  the 
leading  edge  with  a subsonic  region  behind  it.  In  this  case,  the  wing 
leading  edge  is  actually  subsonic  even  though  M = 1.0.  Thus,  for  a 
wing  of  finite  thickness,  it  is  possible  to  have  leading  edge  flow 
separation  for  a supersonic  and  hypersonic  free  stream. 

In  this  investigation,  we  are  concerned  with  the  supersonic  and 
hypersonic  flow  fields  around  a delta  wing  with  supersonic  leading  edges. 
This  study  will  examine  both  the  flow  In  the  boundary  layer  as  well  as 
the  flow  in  the  invlscld  regions.  The  invlscld  flow  field  Is  approxi- 
mately conical  (Refs  13-18)  on  both  sides  of  the  delta  wing.  The 
general  features  of  the  flow  are  shown  In  Figure  2.  On  the  lower  sur- 
face, the  attached  shock  wave  at  the  leading  edge  produces  a cross-flow 


away  from  the  plane  of  symmetry.  Inside  the  Mach  cone,  the  flow  ex- 
pands and  eventually  becomes  parallel  to  the  plane  of  symmetry.  On  the 
expansion  side,  the  region  between  the  Mach  cone  and  the  leading  edge 


contains  attached  lsentroplc  flow.  The  Prandt l-Meyer  expansion  of  the 
free  stream  velocity  component  normal  to  the  leading  edge  causes  a de- 
flection of  the  total  velocity  vector  toward  the  wing  root  chord.  An 
Internal  shock  wave  along  a radial  from  the  vertex  realigns  the  flow 
with  the  plane  of  symmetry  so  that  no  flow  crosses  the  center  plane. 

If  the  internal  shock  Is  weak,  tin'  deflected  flow  remains  attached. 
However,  If  the  shock  wave  is  strong,  the  boundary  layer  will  separate 
resulting  In  a rolled-up  vortex  sheet. 

Since  the  Inviscid  flow  around  the  delta  wing  is  conical,  the  flow 
field  In  any  plane  normal  to  the  wing  chord  has  a similar  solution. 

Thus,  the  entire  flow  field  can  be  represented  by  selecting  a cross-flow 
plane  normal  to  the  axis  of  symmetry.  By  neglecting  the  viscous  effects, 
a qualitative  sketch  of  this  cross-flow  plane  is  shown  in  Figure  3 
(Ref  15). 

On  the  expansion  side  of  the  wing,  the  region  bounded  by  A BOG  is 
hyperbolic.  Elliptical  equations  describe  the  flow  in  the  region  OCBOG. 
The  line  Bll  is  the  cross-flow  sonic  line  which  merges  into  the  internal 
oblique  shock  DKG.  The  Prandt l-Meyer  expansion  is  confined  to  the  area 
ABl)  where  the  line  AB  represents  the  first  ray  of  the  expansion  fan. 

The  line  BEE  is  the  reflection  of  this  expansion  ray  from  t lie  surface 
of  the  Mach  cone  BG  which  emanates  from  the  vertex  of  the  wing.  Ex- 
pansion waves  travel  along  the  straight  characteristics  in  the  Prandt 1- 
Meyer  expansion  and  are  reflected  as  compression  waves  at  the  cross- 
flow  sonic  line.  In  the  region  between  the  cross-flow  sonic  line  and 
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BF,  the  expansion  waves  are  along  one  family  of  characteristics  and  the 
compression  waves  are  along  the  other  family.  The  compression  waves 
below  the  characteristic  AID  are  straight  characteristics  which  are 
eventually  absorbed  by  the  internal  shock  wave.  The  line  OCA  is  the 
wing  surface  and  the  line  COK  represents  the  plane  of  symmetry  (Refs 
14-15). 

For  the  compression  side  of  the  delta  wing,  the  line  AHK  is  a bow 
shock  wave.  The  curvature  of  the  shock  wave  at  KH  makes  the  flow  rota- 
tional in  the  central  region  OKHJ.  The  strength  of  the  shock  is  reduced 
between  the  point  H and  the  plane  of  symmetry.  The  flow  in  the  region 
AHJ  (within  the  bow  shock  wave)  is  isentropic  and  uniform  with  the  area 
JHKO  being  subject  to  the  influence  of  the  wing  vertex.  The  Mach  cone 
or  sonic  line  JH  intersects  the  shock  wave  AHK  at  the  point  H where  the 
shock  begins  to  bend. 

The  spherical  cross-flow  streamline  pattern  for  the  inviscid  flow 
field  is  shown  in  Figure  4.  Outside  the  influence  of  the  delta  wing, 
the  velocity  and  entropy  are  constant  and  the  cross-flow  streamlines 
are  a family  of  straight  lines  that  intersect  at  a point  on  the  plane 
of  symmetry.  However,  when  the  cross-flow  streamlines  cross  a shock 
wave,  these  streamlines  become  curved  lines  that  terminate  at  one  or  more 
singularity  points  called  vortical  singularities.  Those  points  are  cross- 
flow  stagnation  points  where  the  entropy  is  multivalued  (Refs  19-21). 

Thus,  Figures  3 and  4 depict  areas  of  interest  for  which  the  numerical 
solution  of  the  problem  is  desired. 

In  this  description  of  the  supersonic  flow  field,  the  viscous 
effects  and  the  shock  wave-boundary  layer  Interaction  effects  were 
ignored.  The  bow  shock  was  depicted  on  the  compression  side  of  the 
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delta  wing  and  was  attached  at  the  leading  edge.  The  boundary  layer  was 
very  thin  for  high  Reynolds  number  flow  and  thus  the  supersonic  flow 
field  could  be  approximated  as  an  inviscid  conical  flow  field.  However, 
it  was  found  by  experiment  (Refs  14-15)  that  a weak  shock  wave  does 


exist  on  the  expansion  side  of  the  delta  wing  near  the  leading  edge. 
This  very  weak  shock  occurs  because  the  boundary  layer  as  well  as  a 


finite  leading  edge  thickness  prevents  the  compression  side  shock  wave 
from  attaching  itself  at  the  leading  edge,  thus  forcing  it  to  partially 
encircle  the  expansion  flow  field.  A qualitative  sketch  of  this  cross- 


flow pattern  with  viscous  effects  at  low  Reynolds  number  is  shown  in 
Figure  5. 


Hypersonic  Flow  Field.  As  the  flow  is  accelerated  to  hypersonic 
speeds,  the  shock  waves  become  stronger  and  the  boundary  layer  becomes 
thicker.  The  interaction  of  the  inner  shock  with  the  attached  boundary 
layer  is  responsible  for  two  of  the  most  important  features  of  the  lee- 
side  flow  field.  These  arc  the  generation  of  embedded  vortices  in  the 
boundary  layer  at  small  angles  of  attack  (Ref  22)  and  the  shock  induced 
flow  separation  at  large  angles  of  attack  (Refs  17,  18,  23-24).  The 
embedded  vortices  lie  fully  within  the  attached  boundary  layer  and  are 
convected  downstream  through  the  boundary  layer.  The  interaction  of 
these  vortices  with  the  upper  wing  surface  results  in  increased  shear 
and  heat  transfer.  Experimental  measurements  (Refs  13  and  20)  indicate 
that  these  vortices  occur  at  small  values  of  the  viscous  interaction 


parameter 
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For  shock  Induced  separation,  the  shear  layer  rolls  up  into  a pair 
of  conically  growing  vortices  as  seen  in  Figure  6,  These  vortices  re- 
main fully  submerged  in  the  boundary  layer  at  low  Reynolds  number  and 
form  a large  viscous  bubble  inboard  of  the  internal  shock.  With  an 
increase  in  angle  of  attack,  the  internal  shock  and  boundary  layer 
separation  line  move  outboard  until  eventually  the  separation  line 
reaches  the  leading  edge.  The  flow  pattern  in  this  case  resembles  the 
classical  subsonic  leading  edge  separation  pattern  as  shown  in  Figure  1. 
Increased  pressure  and  heat  transfer  rates  are  produced  on  the  leeward 
side  where  the  vortices  impinge  on  the  surface.  Under  some  conditions, 
as  the  embedded  shock  increases  In  strength  the  reattached  flow  under- 
neath the  primary  vortex  might  give  rise  to  a secondary  vortex  as  it 
moves  outward  from  the  centerline  and  separates.  Experimental  measure- 
ments made  by  Cross  (Ref  17)  and  Naravan  (Ref  25)  at  M - 0(1)  exhibit 
the  characteristics  of  shock  induced  separated  flow.  Further  discus- 
sions on  the  flow  field  details  are  presented  in  the  following  sections. 

Experimental  Flow  Field  Studies 

Because  of  the  complexity  of  this  problem,  particularly  at  hyper- 
sonic speeds,  a large  amount  of  experimental  data  exists  on  delta 
wings.  However,  most  of  this  information  is  surface  pressure  .and  heat 
transfer  measurements  (Refs  25-55).  The  majority  of  the  boundary  layer 
and  flow  field  data  comes  from  shadowgraphs,  schlieren  photographs , oil 
flow  pictures,  and  vapor  screen  techniques.  Only  a small  amount  of 
published  detailed  information  has  been  found  for  the  supersonic  and 
hypersonic  flow  over  delta  wings  with  supersonic  leading  edges.  Some 
of  this  detailed  data  consists  of  pitot  pressure  measurements  made  in 
the  plane  of  symmetry  (Refs  22  and  26).  Most  o f this  informat  Ion  consists 
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of  Impact  pressure  surveys  made  in  the  expansion  flow  of  a delta  wing. 
These  pressure  measurements  were  made  by  Fowell  (Ref  13),  Bannink  et  al 
(Refs  14  and  15),  Spurlin  (Ref  16),  Cross  (Ref  17)  and  Monnerie  and 
Werle  (Ref  18).  Fowell  and  Bannink  examined  the  supersonic  inviscid  flow 
field  over  a delta  wing  while  Spurlin,  Cross,  Monnerie  and  Werle  in- 
vestigated the  hypersonic  flow  regime. 

Fowell  (Ref  13)  used  a pressure  rake  to  measure  the  spanwise  pitot 
pressure  above  a flat  delta  wing.  His  measurements  were  made  at  a Mach 
number  of  2.5  and  at  an  angle  of  attack  of  10°  s'.  He  found  that  an  in- 
ternal shock  does  exist  on  the  expansion  side  of  the  delta  wing  and  that 
it  increases  in  strength  as  the  angle  of  attack  is  increased.  The 
strength  of  the  internal  shock  was  quite  small,  being  of  the  same  order 
as  the  weak  bow  shock  on  the  expansion  side  of  the  wing. 

Bannink,  Nebbeling,  and  Reyn  (Ref  14)  did  a more  detailed  investiga- 
tion of  the  supersonic  flow  over  a delta  wing.  Spanwise  pitot  pressure 
measurements  were  made  in  several  planes  perpendicular  to  the  wing  sur- 
face at  various  distances  from  the  wing  vertex.  These  measurements 
focused  on  determining  the  location  and  shape  of  the  internal  shock. 

It  was  found  that  the  internal  shock  was  approximately  conical  and  that, 
away  from  the  wing  surface,  it  bends  towards  the  plane  of  symmetry  of  the 
wing . 

In  1971,  Bannink  and  Nebbeling  (Ref  15)  repeated  their  investigation 
of  supersonic,  inviscid  flow  field  over  a flat  delta  wing.  Tests  were 
conducted  at  a Mach  number  of  2.94  and  at  a Reynolds  number  of  2.65  x 
10^.  The  main  purpose  of  this  study  was  to  determine  the  shape  and 
location  of  the  internal  shock  and  of  the  conical-sonic  line,  the  latter 
being  tiie  locus  of  points  where  the  conical  Mach  number  is  equal  to  one. 
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A comparison  of  pressure  measurements  was  made  in  different  planes 
normal  to  the  root  chord  to  show  that  the  invlscid  flow  is  approxi- 
mately conical.  The  conical  streamlines  were  traced  and  were  found  to 
end  at  the  wing  root  chord.  Thus,  the  conical  stagnation  point  or 
vortical  singularity  is  very  near  the  origin.  The  results  of  this 
investigation  suggest  that  the  Internal  shock  wave  ends  with  zero 
strength  at  the  point  where  the  conical-sonic  line  begins. 

Spurlin,  Cross,  Monnerie,  and  Werle  examined  the  expansion  flow 
over  a sharp  leading  edge  delta  wing  at  hypersonic  speeds.  Monnerie 
and  Werle  (Ref  18)  made  static  and  pitot  pressure  measurements  of  the 
flow  field  at  a Mach  number  of  7.0  and  a Reynolds  number  of  4.5  x 10^. 

The  angle  of  attack  was  varied  from  0°  to  201"  with  flow  separation 
occurring  at  a = 5°.  The  impact  pressure  data  was  presented  as  contour 
curves  on  the  conical  plane.  No  specific  experimental  point  data  was 
given  in  the  reference.  Spurlin  (Ref  16)  made  impact  pressure  measure- 
ments at  a Mach  number  of  10.13  and  a Reynolds  number  of  4.28  x 10". 

This  data  was  used  to  supplement  the  experimental  data  Cross  obtained  for 
the  same  flow  conditions  (Ref  17).  Qualitative  flow  information  was  ob- 
tained from  shadowgraphs,  schlierens,  vapor  screen,  and  oil  flow  pictures. 
Surface  pressures,  surface  temperatures,  and  impact  pressures  were 
measured  at  various  points  in  the  leeside  flow  field  to  determine  the 
magnitude  and  extent  of  the  flow  phenomena.  The  character  of  the 
hypersonic  flow  was  separated  into  three  different  regimes  by  several 
specific  and  dominant  flow  properties.  The  first  regime,  which  occurred 
for  angle  of  attack  in  the  range  cx  < 7°,  was  dominated  by  an  attached, 
quasi-parallel  flow  across  the  surface.  The  second  regime  was  dominated 
by  a conical,  shock-induced  separation  of  the  viscous  layer  for  the  range 
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7°  £ a < 20°.  Finally,  the  third  region  was  identified  by  a base  pres- 
sure induced  separation  of  the  upper  surface  flow  for  a 20°,  which 
eventually  resulted  in  the  separation  of  the  entire  surface  flow. 

Because  of  the  conical  nature  of  the  second  flow  regime,  most  of  the 
comparison  of  numerical  and  experimental  hypersonic  results  was  done 
for  this  regime. 

For  angles  of  attack  from  a = 9°  to  20°,  the  upper  surface  flow 
field  was  dominated  by  a shock-boundary  layer  interaction  that  caused 
a separation  of  most  of  the  wing  boundary  layer.  Separation  of  the 
boundary'  layer  was  caused  by  the  internal  shock  which  was  located  near 
the  Mach  cone.  A viscous  bubble  occurred  near  the  plane  of  symmetry' 
as  a result  of  this  shock-induced  separation.  The  surface' fLow  between 
the  wing  leading  edge  and  the  separation  line  was  attached  and  character- 
ized by  a two-dimensional  Prandtl-Meyer  expansion.  Results  of  the  impact 
pressure  survey  indicated  that  to  a good  approximat ion,  the  flow  was 
conical,  thus  making  possible  a theoretical  or  numerical  analysis  of 
the  flow  in  two  dimensions. 

Theoretical  Flow  Field  Studies 

Several  analytical  techniques  have  been  developed  to  describe  the 
supersonic  and  hypersonic  flow  around  delta  wings.  Most  of  these 
methods  provided  satisfactory  solutions  for  a limited  number  of  flow 
conditions.  Very  few  techniques  have  been  developed  which  can  solve 
the  complete  supersonic  and  hypersonic  flow  field  around  delta  wings 
with  supersonic  leading  edges. 

Supersonic  F 1 ow  Field.  Initial  efforts  at  solving  the  supersonic 
problem  were  directed  towards  using  linearized  conical  flow  theory. 
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Busemann  (Ref  56)  Introduced  this  procedure  and  Stewart  (Ref  57)  and 
Lagerstrom  (Ref  58)  extended  it  to  thin  delta  wings.  This  method  was 
applicable  to  thin  conical  delta  wings  with  weak  shock  waves.  Clarke 
and  Wallace  (Ref  59)  extended  the  linear  theory  to  second  order  in 
incidence,  still  based  on  potential  approximations.  However,  with 
increasing  Mach  number  and  angle  of  attack,  nonlinear  effects  became 
important  and  the  results  of  linearized  theory  were  of  little  value. 

In  fact,  at  high  Mach  numbers  even  the  second  order  theory  was  not 
adequate.  Therefore,  in  order  to  obtain  meaningful  results,  nonlinear 

effects  were  included  in  the  theoretical  analyses. 

» ‘ > , 

One  of  the  first  nonlinear  solutions  to  the  delta  wing  problem  was 
obtained  by  Maslen  (Ref  60) . In  his  treatment  of  the  problem,  he 
assumed  that  the  expansion  side  had  no  shocks  and  thus  was  irrotational . 
Maslen  computed  the  flow  field  characteristics  in  the  hyperbolic  region 
by  using  the  Prandtl-Meyer  relations  and  the  method  of  characteristics. 
The  remaining  subsonic  domain  was  solved  by  using  relaxation  methods. 

On  the  compression  side  of  the  delta  wing,  the  supersonic  flow  properties 
were  calculated  by  using  shock  tables  while  the  remaining  elliptical 
region  was  solved  by  relaxation  techniques.  The  sonic  line  location 
was  determined  empirically  by  using  the  free  stream  Mach  angle  and 
the  shock  wave  angle  from  a wedge.  The  primary  difficulty  encountered 
in  Maslen' s technique  was  that  the  relaxation  procedure  did  not  con- 
verge at  high  supersonic  speeds. 

Fowell  (Ref  13)  improved  on  Maslen's  treatment  of  the  problem  by 
solving  the  governing  equations  for  steady  Inviscid  supersonic  flow  over 
a plane  delta  wing  with  supersonic  leading  edges.  He  showed  that  two 
solutions  exist  for  the  flows  over  the  expansion  surfaces;  a continuous 
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solution  for  flows  below  a critical  a ngle  of  attack  and  a discontinuous 


solution  for  flows  above  this  angle  of  attack.  An  exact  solution  was 
determined  for  the  flow  over  the  compression  surface  and  for  the  con- 
tinuous flow  over  the  expansion  surface.  The  Prandt 1-Meyer  equations 
were  used  to  calculate  the  flow  properties  outside  the  Mach  cone  while 
the  relaxation  techniques  were  used  to  determine  the  flow  properties 
within  the  vertex  domain.  An  Interactive  procedure  was  needed  to 
locate  the  curved  shock  on  the  compress  ion  side. 

For  the  discontinuous  flow  over  the  leeward  side  of  the  delta  wing, 
an  approximate  solution  was  developed  by  Kowell.  He  assumed  the  ex- 
istence of  a planar  shock  normal  to  the  wing  surface  and  of  such  strength 
to  change  the  direction  of  the  outboard  velocity  vector  to  parallel  the 
wing  root  chord.  Two-dimensional  shock  relations  were  used  to  determine 
the  flow  states  on  either  side  of  the  shock  wave.  Kowell  verified  his 
prediction  of  an  internal  shock  wave  by  experimental  study  on  the  ex- 
pans  Ion  side  of  a supersonic  delta  wing. 

It  was  argued  by  Reyn  (Refs  61-63)  that  if  an  irrotational  solution 
for  the  flow  on  the  expansion  side  of  the  delta  wing  did  exist,  this 
wave  pattern  would  lead  to  limit-line  singularities.  This  conical  limit- 
line  would  consist  of  the  line  BFK  (Fig  3)  extending  as  a straight  line 
to  the  wing  surface.  Since  limit-lines  do  not  exist  in  real  flows, 

Bulakh  (Ref  64)  expected  a shock  wave  to  occur  in  the  flow  field.  Bulakh 
introduced  a shock  wave  upstream  of  BKEG  (Kig  3)  and  indicated  that  a 
shock  wave  upstream  of  BC  (Kig  3)  may  also  be  Introduced.  The  numerical 
calculations  following  such  a procedure  were  made  by  Babaev  (Ref  6S) 
and  the  shock  wave  vanished  along  BC  (Kig  3)  and  along  the  characteristic 
BKE  (Kig  3).  A straight  shock  wave  perpendicular  to  the  wing  surface 
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and  extending  to  the  characteristic  BFE  (Fig  3)  was  found.  Babaev  showed 
that  It  was  impossible  to  have  a continuous  flow  solution  over  the  ex- 
pansion side  of  the  delta  wing,  as  described  by  Fowell  (Ref  13). 

In  his  theoretical  development,  Babaev  used  the  1’ rand 1 1 -Meyer  re- 
lations and  the  method  of  characteristics  to  determine  the  hyperbolic 
flow  field  properties  on  the  expansion  side  of  the  delta  wing  (Kef  6 3). 

He  obtained  the  elliptical  flow  field  properties  from  the  irrotational 
potential  flow  theory  bv  using  the  method  of  successive  approximations. 

On  the  compression  side,  Babaev  (Ref  66)  applied  the  oblique  shock  re- 
lations to  the  uniform  flow  bounded  by  the  plane  shock,  the  Mach  cone, 
and  the  wing.  The  flow  properties  in  the  central  region  OKlId  (Fig  3) 
were  calculated  by  the  method  of  successive  corrections  similar  to  the 
method  used  for  the  upper  wing  surface. 

South  and  Klunker  (Ref  67)  improved  on  the  compression  side  flow 
field  results  of  Babaev,  by  using  the  method  of  lines.  In  their  work, 
the  flow  region  was  transformed  into  a rectangular  domain  with  an 
orthogonal  coordinate  grid.  The  shock  wave  and  body  surface  were 
mapped  as  parallel  boundary  lines  in  the  transformed  plane.  On  each 
grid  line,  the  system  of  equations  was  reduced  to  ordinary  differential 
equations  which  were  Integrated  by  a fourth-order  Runge-Kutta  method. 

The  results  of  tills  technique  were  compared  witli  the  numerical  results 
of  Voskresenskl l (Kef  88)  and  were  found  to  be  In  good  agreement. 

Hypersonic  Flow  Field.  For  hypersonic  flow,  most  of  the  theoretical 
analysis  was  focused  on  the  compression  side  flow  field.  Only  two  ref- 
erences were  found  which  attempted  to  calculate  the  leeside  hypersonic 
flow  properties  (Kefs  17  and  68).  In  Cross'  experimental  Investigation 
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(Ref  17),  a simplified  analytical  model  was  developed  to  describe  the 
significant  flow  properties  over  a delta  wing.  Beeman  and  Powers  (Ref 
68)  used  a three-dimensional  method  of  characteristics  to  describe  both 
the  expansion  and  compression  side  flow  fields.  Both  methods  provided 
a fair  agreement  with  experimental  results,  even  though  neither  tech- 
nique included  viscous  and  rotational  effects.  More  significant 
variations  between  theoretical  and  experimental  results  were  noted  in 
regions  of  viscous-invlscid  interaction. 

For  the  compression  side  flow  field,  several  methods  were  used  to 
calculate  flow  field  characteristics.  Thin  shock  layer  theory  was  used 
by  Messiter  (Ref  69)  and  others  to  determine  surface  pressure  distribu- 
tions and  shock  shapes.  This  theory  assumed  that  the  shock  wave  lies 
very  close  to  the  body  and  that  the  pressure  on  the  body  is  determined 
primarily  by  the  shock  shape.  Under  these  conditions,  the  flow  variables 
in  the  intervening  shock  layer  can  be  written  as  a series  expansion  in 
£,  where  £ is  the  density  ratio  across  a plane  shock  lying  at  the  same 
angle  of  attack  as  the  wing.  The  zeroth-order  solution  of  these  series 
gives  the  well-known  Newtonian  results,  modified  possibly  by  the  Busemann 
correction  for  body  curvature.  Messiter  derived  a set  of  approximate 
nonlinear  equations  for  the  conical  flow  and  produced  limited  results 
for  the  lifting  flat  delta  wing  with  a detached  shock.  He  pointed  out 
that  in  the  case  when  the  shock  wave  is  attached  to  leading  edges, 
difficulties  may  arise  in  matching  the  uniform  flow  near  the  leading 
edges  with  the  nonuniform  flow  in  the  central  or  wing  root  region. 

This  work  was  expanded  by  Hida  (Ref  70),  who  made  an  approximate 
allowance  for  wing  thickness.  Squire  (Refs  71  and  72)  and  Shanbhag 
(Ref  73)  obtained  the  exact  numerical  solutions  to  the  governing 
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equations  for  the  flow,  with  detached  shocks,  on  the  compression  surface. 
Their  results  showed  excellent  agreement  with  experiment,  as  did  also 
the  work  by  Htlller  (Refs  74  and  75),  who  extended  the  analysis  to 
include  the  effects  of  yaw.  Squire  (Ref  76),  Woods  (Refs  77  and  78), 

Roe  (Refs  79  and  80),  and  Conor  (Refs  81  and  82)  all  considered  conical 
flows  for  attached  shocks,  again  showing  good  agreement  with  the 
limited  experimental  results  and  exact  solutions  available. 

For  hypersonic  flow  at  low  angles  of  attack  on  the  compression 
side,  the  thin  shock  layer  theory  is  not  valid.  In  this  regime,  linear- 
ized hypersonic  small  disturbance  theory  was  used  by  Malmuth  (Ref  83), 
Ter-Minasyants  (Ref  84),  and  Hui  (Refs  85  and  86)  to  treat  planar  delta 
wings  of  moderate  aspect  ratio.  Malmuth  calculated  the  pressures  on  the 
compression  side  of  a planar  delta  wing  at  infinite  Mach  number.  He 
assumed  that  the  flow  Inside  of  the  Mach  cone  differed  slightly  from  the 
corresponding  two-dimensional  flow  over  a flat  plate  at  the  same  angle 
of  attack.  Based  on  this  assumption,  a linear  perturbation  method  was 
used  to  calculate  the  nonlinear  flow  inside  the  Mach  cone.  Ter- 
Mlnasyants  (Ref  84)  and  Hui  (Refs  85  and  86)  calculated  the  flow  field 
on  both  sides  of  t lie  Mach  cone  at  Mach  numbers  less  than  infinity. 

They  used  a linearized  perturbation  technique  to  calculate  the  flow 
within  the  Mach  cone  and  then  used  l.ighthill's  strained  coordinate 
technique  to  match  this  solution  with  the  exact  solution  of  the  uniform 
flow  outside  of  tlu'  Mach  cone.  When  compared  with  available  exact 
numerical  solutions,  tills  method  gave  almost  identical  results,  except 
near  the  cross-flow  sonic  line,  where  the  numerical  methods  failed  to 
describe  the  discontinuous  flow.  At  the  sonic  line,  this  technique 
predicted  the  discontinuity  and  showed  that  the  discontinuous  pressure 
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slope  has  a square  root  singularity  similar  to  that  described  by  super- 
sonic linear  theory. 

Numerical  Flow  Field  Studies 

Several  references  (Refs  87-93)  were  found  which  used  finite- 
difference  techniques  to  solve  the  supersonic  flow  around  delta  wings. 
Most  of  these  references  presented  a detailed  analysis  and  description 
of  this  conical  flow  problem.  These  references  included  several  papers 
by  Voskresenskii  (Refs  87  and  88),  Kutler  (Refs  89-91)  and  Bazzhin  (Ref 
92). 

Voskresenskii  used  a three-dimensional  implicit  finite-difference 
scheme  to  determine  the  bow  shock  shape  and  the  pressure  distribution 
on  the  lower  wing  surface.  On  the  leeward  side,  the  Prandtl-Meyer 
relations  were  used  to  calculate  the  flow  near  the  leading  edge  while 
the  implicit  numerical  scheme  was  used  to  calculate  the  flow  within  the 
Mach  cone.  In  this  well-posed  initial  value  problem,  Voskresenskii 
approximated  the  initial  shock  wave  shape  and  flow  conditions  near  the 
wing  vertex  and  then  numerically  integrated  the  governing  equations 
until  conical  similarity  conditions  were  sufficiently  satisfied.  This 
technique  produced  a satisfactory  solution  for  both  the  expansion  flow 
field  (Ref  87)  and  the  compression  flow  field  (Ref  8S) , even  though  it 
required  a large  amount  of  computer  storage  and  time. 

Kutler  numerically  computed  the  supersonic,  inviscid,  expansion 
side  flow  over  a planar  delta  wing.  He  used  the  MacCormack  finite- 
difference  scheme  and  a shock  capturing  technique  to  calculate  this 
flow  field.  The  numerical  integration  was  done  in  t he  conical  co- 
ordinate direction,  since  the  governing  equations  were  hyperbolic  in 
that  direction.  His  numerical  results  showed  excellent  agreement  witli 
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the  experimental  data  by  Bannink  and  Nebbellng  (Ref  15)  as  well  as  with 
the  method  of  characteristics  results  by  Beeman  and  Powers  (Ref  68). 

Bazzhin  (Ref  92)  used  the  MacCoraack  finite-difference  scheme  to 
solve  the  time-dependent  governing  equations  for  the  supersonic  flow 
around  a thin  delta  wing  at  low  angles  of  attack.  He  used  a coarse  grid 
in  the  computational  plane  and  thus  was  unable  to  obtain  a true  solution 
in  the  vicinity  of  the  leading  edge.  This  loss  of  accuracy  did  not 
prevent  him  from  obtaining  an  accurate  solution  in  other  areas  of  the 
flow  field.  His  numerical  results  compared  quite  favorably  with  the 
analytical  solutions  by  Voskresenskii  (Refs  87-88). 

Present  Approach 

Although  quite  a few  solutions  were  found  for  the  supersonic  and 
hypersonic  flow  over  delta  wings,  none  of  these  solutions  provided  a 
complete  description  of  the  flow  field.  A number  of  exact  and  approxi- 
mate methods  were  discovered  which  could  accurately  calculate  the 
inviscid  flow  field  around  a delta  wing.  However,  no  satisfactory  theory 
was  found  which  could  describe  the  expansion  side  flow  over  a delta  wing 
with  shock-induced  separation  (Ref  3).  Thus,  the  purpose  of  this  re- 
search was  to  develop  a numerical  solution  for  the  supersonic  and  hyper- 
sonic flow  around  a thin  delta  wing  which  could  adequately  describe  the 
viscous  interaction  phenomena  of  the  flow  with  shock-induced  separation. 

In  order  to  solve  this  problem,  it  was  assumed  that  the  flow  was 
both  conical  and  laminar.  A locally  conical  approximation,  developed 
by  McRae  (Ref  94),  was  successfully  applied  to  the  governing  equations 
for  a compressible,  viscous,  heat  conducting  fluid.  Both  a time- 
dependent  finite-difference  method  and  a shock-capturing  technique  were 
used  to  solve  these  equations.  A computer  program  was  written  to  solve 
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the  governing  equations  In  a generalized  conical  coordinate  plane.  In 
order  to  verify  the  accuracy  of  this  computer  code,  the  supersonic  flow 
around  a cone  was  solved  and  these  results  were  compared  with  both 
experimental  data  and  computed  numerical  results.  The  conical,  viscous 
approximation  was  tested  next  by  computing  both  the  expansion  and  com- 
pression flow  fields  of  a delta  wing  at  supersonic  and  hypersonic  speeds. 
The  results  of  these  calculations,  when  compared  with  experimental  data, 
verified  the  applicability  of  this  technique  for  solving  the  flow  above 
and  below  a planar  delta  wing.  Finally,  the  total  flow  field  around 
a thin  delta  wing  was  solved  at  supersonic  and  hypersonic  speeds  to 
determine  the  interaction  between  the  upper  and  lower  surface  flow 
fields.  These  numerical  calculations  were  compared  with  experimental 
data.  The  results  of  this  research  demonstrated  the  feasibility  of 
obtaining  solutions  to  both  the  supersonic  and  hypersonic  viscous  flow 
fields  over  a planar  delta  wing  by  solving  the  conical  governing 
equations. 
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II.  Governing  Equations 


In  this  chapter,  the  unsteady  governing  equations  for  a viscous, 
compressible  fluid  are  defined  along  with  the  basic  fluid  flow 
assumptions.  The  equations  are  transformed  into  generalized  conical 
coordinates  by  using  the  relations  in  Appendix  A.  A conical  approx- 
imation is  applied  to  the  resulting  equations  and  then  these  equations 
are  put  in  nondimenslonalized  conservative  form. 


Navier-Stokes  Equations 

The  fundamental  governing  equations  for  a compressible,  viscous, 
heat  conducting  fluid  in  Cartesian  tensor  form  are  given  by 
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In  addition  to  the  above  equations  which  express  the  conservation 
of  mass,  momentum,  and  energy,  it  is  necessary  to  include  a state  equa- 
tion which,  for  a perfect  gas,  is  given  by 

f - PRT  (9) 


The  viscosity  coefficient  p Is  a function  of  temperature  and  is  ade- 
quately approximated  by  Sutherland's  serai-empirical  equation 


- k.  (f)*  [ilVM 

h + c0/TJ 


(10) 


where  the  constand  Cq  has  been  experimentally  determined  to  be  198.6  R 

for  air.  The  Sutherland  law  takes  into  account  the  weak  attraction 

field  surrounding  the  molecule  and  therefore  provides  a more  rapid 

viscosity  variation  than  the  rigid  sphere  result  of  kinetic  theory 

\< 

which  gives  the  result  U - T . 

Equations  2 through  10  provide  seven  basic  equations  in  the  seven 
unknowns,  u,  v,  w,  p,  T,  p,  and  p,  if  the  bulk  viscosity  coefficient 
is  set  equal  to  zero.  This  implies  that  the  second  coefficient  of 
viscosity  X is  equal  to  - ^3  p.  For  a monatomic  gas,  this  step  is 
Justified  since  the  bulk  viscosity  coefficient  vanishes  when  the 
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molecule  has  no  Internal  degrees  of  freedom.  However,  for  a poly- 
atomic gas,  the  bulk  viscosity  effect  is  not  always  zero  and  can  be  of 
the  same  order  of  magnitude  as  the  molecular  viscosity  in  sound  propa- 
gation and  shock  wave  structure.  A more  detailed  discussion  of  the 
significance  of  bulk  viscosity  is  given  by  Vincenti  and  Kruger  (Ref  95) 
In  the  present  investigation  it  is  neglected  since  the  numerical  grid 
spacing  used  in  the  shock  transition  region  is  too  coarse  to  allow 
resolution  of  the  shock  structure.  This  lack  of  resolution  causes 
numerical  difficulties  which  will  be  discussed  in  chapter  three. 

There  are  several  approaches  which  may  be  used  to  derive  the 
governing  equations.  However,  the  most  commonly  used  method  is  the 
continuum  approach.  In  this  approach,  a volume  of  arbitrary  dimensions 
is  fixed  in  space  and  enclosed  by  an  imaginary  surface  through  which  a 
fluid  is  flowing.  The  continuity  equation  is  derived  from  the  law  of 
conservation  of  mass  which  states  that  mass  can  neither  be  created  or 
destroyed  (Eulerian  Method).  Therefore,  the  time  rate  of  change  of 
mass  within  the  volume  must  equal  the  net  inward  rate  at  which  mass 
crosses  the  surface.  The  Navier-Stokes  or  momentum  equations  are 
applications  of  Newton's  law  of  motion  which  equates  the  total  force 
acting  on  the  mass  in  the  arbitrary  volume  to  the  rate  of  change  of 
linear  momentum.  The  total  force  consists  of  the  pressure,  the  body 
force,  and  the  viscous  stress  which  is  assumed  to  be  linearly  related 
to  the  rate  of  strain.  The  energy  equation  is  based  on  the  first  law 
of  thermodynamics  which  states  that  the  increase  in  internal  energy  of 
the  mass  of  fluid  in  the  volume  is  equal  to  the  sum  of  the  heat  added 
and  work  done  on  the  fluid  mass.  In  the  absence  of  chemical  reaction 
and  radiation  heat  flux,  the  heat  added  is  due  to  conduction  and  is 


given  by  the  Fourier  law  of  heat  conduction.  The  work  done  on  the  fluid 


Is  that  due  to  the  pressure  and  shear  stress  acting  on  the  volume  sur- 
face. A more  detailed  discussion  of  this  approach  of  deriving  the 
governing  equations  may  he  found  in  Schlichtlng  (Ref  96),  Yuan  (Ref  97), 
and  Chow  (Ref  98). 

In  this  Investigation,  the  following  assumptions  are  applied  to 
the  governing  equations: 

1.  The  specific  heats,  c and  c , are  both  constant  and  thus  the 
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specific  heat  ratio  y =*  ^ Is  also  constant. 
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The  Prandtl  number  Pr  » u c /k  is  constant  throughout  the 
a given  value  not  necessarily  unity. 

The  coefficient  of  thermal  conductivity  is  computed  from 
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4.  Dissociation  and  ionization  effects  are  not  present  in  the 


flow. 


5.  Stokes'  relationship  between  the  first  and  second  viscosity 
coefficients  is  valid. 

6.  There  is  no  external  heat  addition. 

7.  The  body  forces  are  negligible. 

8.  There  are  no  chemical  reactions  and  radiation  heat  flux. 

From  the  above  assumptions,  it  can  bo  seen  that  the  real  gas 

effects  associated  with  high  Mach  number  and  high  temperature  flows 
arc  neglected.  However,  for  comparison  with  wind  tunnel  data,  the 
calculated  impact  pressures  are  corrected  for  real  gas  effects  by  using 
the  caloric  imperfection  charts  in  NACA  Report  1135. 
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Conservative  Form 


In  order  to  process  the  discontinuities  in  the  flow  field  correctly, 
the  conservative  form  of  the  governing  equations  is  required.  This  con- 
servative form  was  initially  used  by  Courant  and  Friedrichs  (Ref  99)  for 
inviscid,  compressible  flow.  Later  Lax  (Ref  100)  applied  this  form  to 
a finite-difference  scheme.  Gary  (Ref  101)  and  Abbett  (Ref  102)  showed 
that  using  the  non-conservative  equations  produced  significant  errors 
in  shock  speed  and  location,  while  Longley  (Ref  103)  showed  that  the 
conservative  form  yielded  the  correct  shock  speed  for  a variety  of 
finite-difference  schemes. 

When  the  governing  differential  equations  are  placed  in  conserva- 
tive form,  the  independent  variables  become  p,  pu,  pv,  pw,  and  pe . The 
use  of  these  variables  in  conservative  finite-difference  schemes 
assures  the  conservation  of  mass,  momentum,  and  energy  across  a shock 
wave.  This  can  be  easily  demonstrated  by  considering  a stationary 
normal  shock  wave.  The  truncation  error  of  a finite-difference  scheme 
depends  on  the  size  of  the  higher  derivatives  in  the  Taylor  series 
expansions  for  the  differentials.  In  the  variables  p,  u,  v,  w,  and  T, 

I 

the  shock  causes  discontinuities  in  the  continuum  solution,  but  the 
solution  is  continuous  in  several  of  the  conservative  variables.  Since 
Rankine-Hugoniot  relations  for  normal  shocks  are  based  on  global  con- 
servation of  fluid  properties,  these  equations  are  independent  of  the 
details  within  the  shock  structure.  As  a result,  the  conservative 
finite-difference  techniques  applied  to  the  governing  equations 
satisfy  the  Rankine-Hugoniot  relations  and  thus  predict  the  correct 
shock  conditions  across  a shock. 

The  continuity  equation  is  already  in  conservative  form.  The  g 
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Con  t o Coord  i nat  os 

The  general  transformation  from  tlio  physical  plane  (x,  y,  z)  to 
the  conical  transformed  plane  (i.,  r),  f.)  is  given  by 
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Assuming  that  an  Inverse  exist,  then 
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Thus,  the  Jacobian  determinant  is  given  by 
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It  is  important  that  the  transformation  equations  define  a mapping 
which  establishes  a one-to-one  correspondence  between  each  point  in  the 
physical  plane  to  one  and  only  one  point  in  the  transformed  plane,  and 
conversely.  This  condition  is  satisfied  if  (a)  the  transformation 
functions  are  continuously  differentiable  and  (b)  if  the  Jacobian  of  the 
transformation  is  nonsingular  at  each  point.  When  these  conditions 
occur,  there  is  a sphere  Nq  about  the  point  such  that  the  inverse  func- 
tion exists  for  all  x,  y,  z in  Nq.  Thus,  the  transformation  is  guaran- 
teed in  a local  fashion  only,  and  the  functions  in  Eq  16  must  be  select- 
ed so  that  they  possess  desirable  characteristics  in  the  region  of 
interest . 

By  applying  the  chain  rule  of  differential  calculus,  the  differen- 
tial operators  become 
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A comprehensive  set  of  transformed  derivatives  is  given  in  Appendix  A. 

Once  the  transformed  coordinates  are  generated  for  a given  physical 
domain,  the  set  of  partial  differential  equations  and  their  associated 
boundary  conditions  are  transformed  utilizing  the  relations  given  in 
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Appendix  A.  Substitution  of  the  above  differential  operators  into  the 
conservative  Cartesian  equations  transforms  the  governing  equations  into 
the  (4,0.4)  coordinate  system.  The  resulting  equations  in  weak  conser- 
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Conical  Approx imat ion 

The  definition  for  supersonic  conical  flow,  given  by  Liepmann  and 
Roshko  (Ref  106),  is  one  in  which  the  fluid  properties  are  invariant 
along  each  ray  emanating  from  a vertex.  An  essential  element  which 
follows  from  this  definition  is  that  conical  flow  problems  do  not  in- 
volve a characteristic  length.  A body  is  called  conical  if  its  surface 
is  made  up  of  rays  from  a vertex.  Boundary  conditions  are  conical  if 
the  boundary  values  of  the  fluid  properties  are  constant  along  rays 
from  an  origin.  All  natural  features  of  the  flow  field  are  conical 
such  as  embedded  shocks,  bow  shocks,  and  contact  surfaces. 

The  basic  theorem  of  conical  flow  is  "a  flow  field  is  conical  if 
the  boundary  conditions  are  conical."  This  fundamental  theorem  applies 
to  linearized  potential  flows  as  well  as  compressible,  invisetd  flows. 
It  does  not  hold,  however,  for  viscous  fluids.  The  theorem  can  be 
demonstrated  quite  easily  by  applying  dimensional  analysis.  Since  a 
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conical  boundary  does  not  contain  a characteristic  length,  then  it  can 
be  shown  that  no  combination  of  fluid  properties  given  at  infinity  can 


be  constructed  to  give  a dimension  of  length.  This  is  not  true,  if 
viscosity  is  considered  since  [V,/u]  * [ L-],  Thus,  the  boundary  layer 
over  a delta  wing  is  a function  of  the  distance  from  the  vertex. 

In  this  investigation,  a conical  approximation  is  applied  to  the 


governing  equations.  It  is  assumed  that  the  viscous  region  thickness 
is  much  smaller  than  the  length  scale  L of  the  wing.  The  leading  edge 
of  the  delta  wing  is  assumed  to  be  sharp  and  the  boundary  layer  is 
laminar  and  attached  at  the  leading  edge  (see  Fig  7).  Since  the 
viscosity  effect  in  the  boundary  layer  introduces  a length  scale,  which 
appears  only  in  the  Reynolds  number  and  radial  derivatives,  then  at 
large  Reynolds  numbers,  the  gradients  in  the  radial  direction  are  much 
smaller  than  those  in  the  cross  flow  directions  (0,4) • Hence, 
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This  is  the  mathematical  definition  of  a locally  conical  flow.  Grad- 
ients do  exist  in  the  radial  direction,  but  they  are  sufficiently  small 
to  be  neglected.  When  this  requirement  is  imposed  on  the  generalized 
transformed  coordinates,  it  becomes 
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Thus,  the  flow  properties  p,  u,  v,  w,  p,  and  T are  approximately  invar- 
iant along  rays  of  constant  r,(r).  However,  it  should  be  noted  that  the 
Reynolds  number  is  a function  of  5 and  therefore  a characteristic 
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length  must  be  Included  in  the  computation. 

Anderson  (Ref  107)  was  the  first  to  introduce  this  conical  approxi- 
mation and  McRae  (Ref  94)  later  demonstrated  it  in  his  calculations  of 
supersonic,  viscous  flows  over  a cone.  From  experimental  data  by  Cross 
(Ref  18)  and  from  discussions  with  Barber  (Ref  31),  it  was  found  that 
the  viscous  layer  is  laminar  and  approximately  conical  under  certain 
hypersonic  flow  conditions.  For  high  Reynolds  number  supersonic  flow, 
the  boundary  layer  is  so  thin  that  its  effect  on  the  total  flow  field 
is  minimal.  Thus,  a locally  conical  flow  approximation  for  the  viscous- 
inviscid  flow  field  over  a delta  wing  is  adequate  to  model  the  test 
cases  being  considered. 


Nondimensionalization 

Before  numerically  solving  the  resulting  governing  equations,  it 
is  advantageous  to  write  these  equations  and  other  relationships  in 
nondimensional  form.  The  pressure  and  density  are  made  dimensionless 
with  respect  to  free  stream  stagnation  conditions  while  velocity  and 
total  energy  are  nondimensionalized  with  respect  to  the  maximum  adiabatic 
velocity 
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The  normalizing  length  scale  L is  the  characteristic  length  used  in 
determining  the  free  stream  Reynolds  number. 

The  following  nondimensional  relationships  are  used  in  the  govern- 
ing equations: 
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and  where  the  subscript  "o"  refers  to  free  stream  isentropic  total 
condit ions . 

We  also  define  a reference  Reynolds  number  as 
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and  a Prandtl  number  as 
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If  we  drop  the  prime  notation  with  the  understanding  that  all  quantities 
arc  nondimens ional  values  unless  otherwise  indicated,  then  the  final 
form  of  the  governing  equation  is 

0 (29) 
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The  nondimensionalized  heat  flux  and  stress  terms  are 
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The  normalized  equation  of  state  becomes 
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and  the  nondimonsional  form  of  the  Sutherland  viscosity  formula  is 
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Coordinate  Trans format  tons 

Although  this  problem  is  developed  for  generalized  conical  coordin- 
ate transformations,  only  two  specific  coordinate  transformations  are 
used.  These  include  a cylindrical-type  coordinate  transformation  for 
flow  over  a cone  and  a conical  coordinate  transformation  for  flow  over 
a delta  wing. 

For  the  flow  field  calculations  over  a cone,  the  coordinate  trans- 
formation is 


The  physical  and  computational  planes  are  shown  in  Figures  8a  and  8b, 
respectively.  The  transformation  derivatives,  which  are  used  in  the 
calculations,  are  developed  in  Appendix  A. 

For  the  flow  field  around  a thin,  planar  delta  wing,  a constant 
spaced  grid  is  used  in  the  physical  plane.  The  conical  coordinate 
transformation  is 

C-x  r\  = Z £ « - (47) 

X X 

where  the  physical  and  computational  domains  are  depicted  in  Figures 
9a  and  9b,  respectively.  The  transformation  derivatives  are  given  in 

Appendix  A. 

These  coordinate  transformations  and  their  derivatives  are  applied 
to  the  conical  governing  equations  for  a compressible,  laminar,  viscous 
flow.  The  equations  are  solved  by  a numerical  finite-difference  proce- 
dure which  Is  discussed  in  the  next  chapter. 
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III.  Numerical  Methods 


In  this  chapter,  several  basic  concepts  of  numerical  methods  are 
discussed  for  solving  time-dependent,  multi-shock  flow  field  problems. 
The  MacCormack  finite-difference  scheme  is  examined  and  then  applied  to 
the  fundamental  governing  equations.  Several  numerical  models  are  used 
to  represent  the  initial  and  boundary  conditions.  The  requirements  of 
stability,  consistency,  and  convergence  are  derived  and  discussed  for 
solving  the  resulting  numerical  equations. 

Basic  Concepts 

For  this  conical  flow  problem,  the  steady  state  solution  is  sought 
by  numerically  solving  the  unsteady  flow  equations.  This  steady  state 
solution  is  determined  when  all  the  time  derivative  terms  in  the  govern- 
ing equations  vanish  in  the  limit  of  large  times.  Such  a procedure  for 
obtaining  the  steady  state  solution  is  sometimes  referred  to  as  the 
"asymptotic"  method  (Ref  108) . This  method  was  advocated  by  such  early 
proponents  as  J.  von  Neumann  and  R.  Rlchtmver  (Ref  109)  because  it 
represents  the  physical  situation  more  naturally  than  the  method  of 
solving  the  steady  state  governing  equations.  The  successive  iterations 
needed  to  determine  the  "steady"  state  solution  may  be  looked  upon  as 
the  evolution  of  the  flow  field  with  time. 

Under  the  conical  flow  approximation,  the  radial  derivatives  of  the 
flow  properties  can  be  set  equal  to  zero.  This  reduces  the  number  of 
independent  variables  in  the  transformed  plane  to  three,  which  include 
time  and  the  two  cross-flow  spatial  coordinates  ant*  ^ ( x'x  ) 

This  conical  flow  problem  is  then  solved  as  an  initial  value  problem. 
Initial  conditions  are  specified  at  t=0  along  with  appropriate  boundary 
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| conditions,  and  the  solution  is  numerically  integrated  in  time  until 

steady  state  conditions  are  reached.  Elimination  of  £(r)  as  one  of 
the  independent  variables  in  the  problem  reduces  the  amount  of  storage 
space,  numerical  computation,  and  resulting  computer  time  of  the 
problem. 

i 

By  applying  finite-difference  techniques,  the  partial  differential 
operators  in  the  governing  equations  are  replaced  by  finite-difference 
quotients.  This  allows  the  governing  equations  to  be  approximately 
represented  by  a set  of  difference  equations.  Solutions  to  the  dif- 
ference equations  are  obtained  at  the  intersection  points  of  the  com- 
putational grid  lines.  These  nodal  points  may  be  specified  by 

a double  subscript  (i,j)  and  a superscript  n,  where  n refers  to  times 
t=nAt  and  where  At  is  the  time  step  that  the  solution  is  advanced 

during  each  cycle. ' All  the  dependent  variables  at  the  nodal  points  are 

I 

initially  prescribed  as  free  stream  values  and  the  unknown  values  at 
t>0  are  calculated  during  numerical  integration.  Figure  10  illustrates 
the  mesh  of  nodal  points  in  the  transformed  plane. 

There  are  several  types  of  finite-difference  quotients  which  may 
be  used  to  represent  the  partial  derivatives.  An  implicit  difference 
quotient  approximates  the  partial  derivative  at  an  advanced  point  in 
time  and  is  written  so  that  the  dependent  variable  is  expressed  in 
terms  of  its  neighboring  values  and  known  previous  values.  This 
implicit  representation  requires  the  solution  of  a set  of  simultaneous 
equations  in  order  to  calculate  the  dependent  unknown  as  each  time  step 
is  taken.  The  implicit  form  lias  the  advantage  of  universal  stability 
for  t>0,  however  the  method  requires  a more  complicated  computational 
procedure. 
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The  explicit  differential  quotient  represents  a partial  differential 
operator  with  only  one  unknown.  The  explicit  form  is  solved  directly 
(explicitly)  in  terms  of  known  previous  values  and  boundary  conditions. 

As  each  step  in  time  is  taken,  the  unknown  quantities  are  calculated  one 
at  a time.  The  explicit  difference  quotient  has  the  advantage  of  pro- 
gramming simplicity,  although  there  are  stability  restrictions  on  its 
use.  Because  of  Lts  successful  application  to  compressible,  viscous 
multiple  shock  flows,  only  the  explicit  finite-difference  quotients  are 
used  in  this  investigation. 

Two  distinct  approaches  are  used  to  compute  multishock  flow  fields. 
One  is  referred  to  as  a shock-fitted  technique  and  the  other  as  a shock- 
capturing method.  For  shock-fitting,  the  shock  is  transformed  onto  a 
coordinate  surface  and  treated  as  a discontinuity.  The  values  at  the 
shock  surface  are  determined  by  explicitly  solving  the  Rank ine-Hugoniot 
relations.  Although  it  is  conceptually  possible  to  treat  all  shocks  as 
discontinuities,  this  method  involves  complex  logic  to  fit  all  the  in- 
ternal and  bow  shocks. 

The  shock-capturing  technique  is  a method  by  which  multiple  shocks 
can  be  accurately  determined  numerically  without  knowing  initially  where 
the  shocks  will  form.  Shock  waves  are  automatically  allowed  to  form  and 
decay  without  employing  any  special  shock-fitting  procedures  which  in- 
volve explicit  use  of  the  Rankine-llugonlot  equations.  The  shock- 
capturing technique  requires  that  the  governing  equations  be  differenced 
in  weak  conservative  form.  Shocks  which  exist  or  do  form  are  generally 
spread  over  3 or  4 mesh  intervals.  Thus,  a finer  mesh  interval  is  re- 
quired in  order  to  specify  the  shock  location.  This  method  also  has  a 
tendency  to  give  spurious  fluctuations  of  flow  field  quantities  In  the 
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vicinity  of  the  shock.  These  fluctuations  can  be  dampened  by  use  of  an 


artificial  viscosity  term.  Unfortunately,  additional  shock  smearing 
occurs.  However,  the  results  obtained  by  using  this  method  have  been 
shown  to  be  quite  accurate  when  compared  to  shock-fitting  methods 
outside  of  the  immediate  shock  location.  Therefore,  in  this  investiga- 
tion shock-capturing  techniques  are  used  to  compute  the  conical  flow 


field  around  a delta  wing. 

Finite-Difference  Scheme 

The  selection  of  the  numerical  shock-capturing  algorithm  involves 
considering  such  factors  as  accuracy,  computational  speed,  storage  re- 
quirements, and  programming  complexity.  Unfortunately,  there  is  no 
generally  accepted  optimal  scheme,  even  when  the  order  of  the  numerical 
approximation  is  specified.  A survey  of  various  first  and  second-order 
schemes  by  Kutler  (Ref  110)  shows  that  the  highly  dissipative  nature  of 
the  first  order  accurate  methods  makes  them  generally  unsuitable  for 
shock-capturing  calculations.  Among  the  second-order  schemes,  the 
predictor-corrector  combination  oy  MacCormack  was  found  to  be  the  best. 
Anderson  (Ref  111)  compared  the  second-order  MacCormack  method  with  both 
the  third-order  Rusanov  and  Kutler-Lomax-Warming  methods.  He  applied 
these  numerical  methods  to  a set  of  one-dimensional  hyperbolic  equa- 
tions 

|f  + = ° (48) 

3t  3x 

where  F.  and  F are  n-component  vectors  and  x is  an  arbitrary  direction. 

By  utilizing  von  Neumann's  stability  analysis,  the  stability  criterion 
for  these  equations  was  found  to  be 

I X I < 1 (49) 
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where  A is  the  Courant  number  given  by 
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and  O represents  the  maximum  eigenvalue  of  the  Jacobian  matrix  -r-r  . 
Anderson  discovered  that  the  second-order  MacCormack  method  gave 
acceptable  results  for  most  cases.  For  a Courant  number  of  approxi- 
mately one,  the  second-order  MacCormack  scheme  provides  the  best  re- 
sults from  among  the  examples  tested.  However,  when  the  Courant 
number  varies  appreciably,  third-order  schemes  provide  better  results. 

In  certain  critical  flow  regions,  such  as  multiple  shock  flows  In  which 
shock  waves  of  various  strengths  are  located  close  together,  the  third- 
order  schemes  are  capable  of  obtaining  results  where  second-order 
methods  would  fail  (Ref  112).  Since  these  critical  flow  regions  re- 
quiring third-order  methods  are  not  expected  to  be  encountered  in  this 
study,  MacCormack's  method  was  selected  for  the  numerical  Integration  of 
he  flow  field. 

The  second-order  MacCormack  scheme  is  a non-centered , preferen- 
tial, predictor-corrector  algorithm  (Ref  113)  which  has  been  used  for 
a wide  variety  of  flow  calculations.  The  method  is  a variant  of  the 
two-step  Lax-Wendroff  method  (Ref  114),  but  because  of  its  non-centered 
nature,  the  dependent  variables  need  not  be  evaluated  at  half-mesh  in- 
tervals. When  the  MacCormack  two-step  explicit  method  is  applied  to 
the  normalized  conservative  equations 


3U  , 9F  . 3C  , ..  _ 

ae  + + K + " 0 

the  following  predictor-corrector  equations  result: 


(51) 
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Predictor 


U?+J  - U"  - ^ ! (1-e  ) F"  - (l-2e  ) F" 

i.j  i.j  An  | n i+i* j n i,j 


e f.  . , 
n i-i.j 


At 

A£ 


(l-ec)  G 


C'  “i.j+l 


(52) 


( 1— 2e _)  g"  - er  g" 


- At  H 


i.j 


Corrector 


n+l  1 

i.j  ' 2 


or  , + or':  - tt  [ e_  f 


J?  . 

i.J 


n+l  jit 

i.j  An 


n+l 

n *i+i,j 


where 


+ <1-2en>  + (VU  F 


n+l  "I 

i-l.j  J 


- t-[Ec  Cu  + (1-2£s>  °Z 


(53) 


n = iAn  £ = jA£  t = nAt 


(54) 


and 


Fj  . - F(U"  ) T*+1.  = F(uf  b , etc 


(55) 


This  scheme  permits  four  possible  variations  for  replacing  the  time 

derivatives  in  the  predictor  and  corrector  steps: 

I:  e =0  = 0 

H £ 


II:  e_  *=  1 e_  = 1 


III:  e_  = 0 e. 


IV:  e_  = 1 c. 


(56) 
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Each  version  yields  a slightly  different  result.  For  this  reason,  it 
is  termed  a preferential  scheme.  MacCormack  (Ref  113)  suggested  that 
the  four  combinations  be  cyclically  permuted  in  order  to  obtain  the 
most  unbiased  result.  However,  it  has  been  shown  that  for  shocks  which 
propagate  away  from  a body  surface  in  the  computational  piano,  method  I 
tends  to  yield  the  best  results  (Refs  110  and  115) . Since  the  shocks 
in  this  flow  field  problem  originate  at  the  body  surface,  method  I was 
used  for  the  numerical  integrations.  Thus,  the  following  predictor- 
corrector  equations  were  used. 

Predictor 


= u"  . 

- AM 

Fn 

- Fn  ( 

An  | 

1+1,  j 

1 

At 
~ AC 

| Gi,j+1 

1 

_ rn  | 

G±,i  i 

1 -4t 

Corrector 

un+1 

±,i 

1 

Vh  + 

pn+1  At  rFn+l_F: 

i,j  An  |_ 

2 

At 
~ AC 

rG^  _ 

LGij 

cn+1 

i.J-1  J 

- At  Hn+1 

1 , J 

(57) 


(58) 


In  this  version,  forward  differences  are  used  in  the  predictor  step  to 
3F  3G 

approximate  -r—  and  derivatives  while  backward  differences  are 

3n  at, 

used  for  these  terms  in  the  corrector  step.  This  method  minimizes  the 
post  and  precursor  oscillations  near  the  shock. 

In  order  to  evaluate  the  shear  stress  and  heat  flux  terms  appearing 
in  the  F and  C matrices,  backward  differences  are  used  in  the  predictor 
step  and  forward  differences  in  the  corrector  step.  A central  dif- 
ference operator  is  used  for  and  qj  in  the  H matrix.  This  dif- 
ferencing procedure  results  in  a central  difference  approximation  with 
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: 


second-order  accuracy  for  all  shear  stress  and  heat  flux  terms  (Ref  116) 

For  example,  if  we  consider  the  F matrix  containing  the  velocity 
9u 

gradient  gpp  , then  the  predictor  term  is 


— T Fn  - Fn  • 

An  L i+l.J  1J  J 


(59) 


and  the  velocity  gradient  terms  can  be  written  as 


and 


n n 

Ui±l,j  ~ Ui,j 

An 


n n 

Ui,j  ~ ui-l,j 

An 


for  the  F_^+^  ^ term 


for  the  F . . term 


For  the  corrector  term,  we  have 


— r Fn+i  - Fn+i  * 

An  [Fi,j  Fi-l,j  J 


where  the  velocity  gradient  terms  are 


^ the  F"+1  term 

1+1 >J  i,J 


and 


An 


Ul,.1  ~ “i-1.1  f°r  the  F£l,j  term 
An 


(60) 


(61) 


(62) 


(63) 


(64) 


This  results  in  a f inite— dif ference  quotient  of  second— order  accuracy 
_ 92u 

for  the  derivative  centered  at  (i,j). 


When  the  cross  derivative  is  calculated  l.e.,  inside  the  F 

matrix,  then  the  velocity  gradient  terms  for  the  predictor  are 


f 
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This  approximation  results  in  a finite-difference  quotient  of  second- 

32u 

order  accuracy  for  the  derivative  centered  at  (i,j). 

By  using  these  finite-difference  equations,  the  calculations  are 
advanced  in  time  from  the  initial  conditions  until  the  steady  state 
solutions  are  achieved.  At  each  time  step  of  the  integration  process, 
it  is  necessary  to  calculate  the  vectors  F,  G,  and  H,  which  are  func- 
tions of  the  dependent  and  independent  flow  variables.  The  dependent 
normalized  flow  variables  are  obtained  at  each  interior  grid  point  by 
decoding  the  U vector 


U1 

P 

U2 

Pu 

U3 

= 

Pv 

U4 

Pw 

U5 

Pe 

r 


In  this  decoding  procedure,  the  primitive  flow  variables  are  solved  for 


in  the  following  manner 


P“U! 

u - U2/U1 

v - U3/Ux 

w - U4/Ux 

e - U5/U1 

P - 2yp  (e  - 

T - - 
P 


u2  + v2  + w2) 


After  the  decoding,  the  boundary  conditions  are  determined  on  the  body 
surface  and  at  the  plane  of  symmetry.  Thus,  with  the  free  stream  and 
boundary  conditions  as  the  initial  starting  solution,  Eq  51  is  integra- 
ted with  respect  to  time  and  the  flow  variables  are  advanced  from  t=t” 
to  t*tn  + At.  If  the  flow  variables  in  the  t=tn  plane  and  the  t=tn  + At 
plane  do  not  satisfy  the  convergence  criteria,  then  the  computed  flow 
variables  at  t=tn  + At  are  transferred  to  the  t=tn  plane  and  become  the 
new  starting  solution.  This  iteration  process  is  continued  until  there 
is  little  change  between  the  flow  variables  at  t=tn  + At  for  two  sub- 
sequent cycles.  At  this  point,  the  solution  is  assumed  to  have  con- 
verged. The  following  inequality  is  used  as  the  termination  criteria 


,n  ,n+l 

*j.j  - +i.i 


for  l<i<i  and  l<j<j 

max  max. 


where  <J>  . represents  each  primitive 

i * J 


variable  at  t=t  , ^ is  a constant,  and  e is  the  convergence  criteria 

(typically  10  . Further  discussions  on  the  concept  of  iterative 


52 


convergence  and  the  determination  of  iteration  errors  can  be  found  in 
Appendix  B. 


Boundary  Conditions 

An  important  aspect  in  the  successful  application  of  thq  MacCormack 
finite-difference  method  is  the  proper  treatment  of  the  boundary  con- 
ditions. There  are  numerous  schemes  available  for  simulating  the  various 
types  of  physical  boundaries  encountered  in  a fluid  flow  problem  (Refs 
102,  117,  118).  However,  there  is  no  rigorous  mathematical  theory  con- 
cerning which  numerical  boundary  conditions  to  select  in  order  to  get  a 
unique  solution.  Thus,  one  must  examine  the  physical  flow  field 
situation  and  take  into  account  the  mathematical  nature  of  the  governing 
equations  in  order  to  select  which  numerical  boundary  conditions  to 
apply. 

In  this  investigation,  there  are  four  types  of  boundary  conditions 
which  require  numerical  modeling.  These  are: 

1.  Body  Surface 

2.  Symmetry  Plane 

3.  Free  Stream  Boundaries 

4.  Singularity  Points 

The  finite-difference  equations  used  to  describe  these  boundary  con- 
ditions are  developed  in  the  following  subsections. 

Body  Surface  Boundary  Conditions.  On  the  surface  of  the  delta  wing 
and  the  cone,  the  boundary  conditions  are 

u * v * w *>  0 (72) 


and 


T-T 


wall 


or 


9T 


- 0 


(73) 
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Equation  72  guarantees  that  the  velocity  component  tangent  to  the  sur- 
face satisfies  the  viscous  no-slip  conditions  and  that  the  velocity 
component  normal  to  the  body  surface  is  zero.  The  thermal  conditions 
on  the  surface  are  either  isothermal  or  adiabatic,  as  indicated  by  Eq 
73.  The  isothermal  boundary  condition  requires  that  a temperature  be 
specified  at  the  surface  nodal  points.  The  adiabatic  boundary  con- 
dition requires  that  a first  or  second-order,  one-sided,  finite-dif- 
ference term  be  used  to  model  the  zero  temperature  derivative.  In  this 
investigation,  all  surfaces  are  isothermal. 

The  surface  pressure  or  density  cannot  be  obtained  from  the  boundary 
conditions,  and  thus  must  be  calculated.  There  are  several  methods  for 
calculating  the  surface  pressure  or  density;  however,  only  two  methods 
seem  preferable.  Either  the  continuity  equation  must  be  satisfied,  for 
which  an  iteration  for  the  density  at  the  body  surface  is  accomplished, 
or  the  normal  momentum  equation  must  be  satisfied  for  the  component  of 
the  pressure  gradient  normal  to  the  surface. 

The  first  technique  consists  of  calculating  the  density  at  the 
surface  from  the  continuity  equation 


9p  3E  9u  3v  3w 

3t  wall  -pai  as  -p3y  as  ^ as 


(74) 


where  the  velocity  and  the  gradients  parallel  to  the  surface  are  zero. 

This  method  has  the  advantage  of  being  simple  and  convenient  when  using 

suitable  second-order  one-sided  difference  quotients  for  the  spatial 

derivatives  and  a first-order  forward  time  difference  quotient  for 

■1^-  . However,  this  technique  may  lead  to  divergence  of  the  numerical 
ot 

results  for  t-*°°.  This  is  particularly  true  for  separated  flows  and 
flows  near  a leading  edge. 


L 


: 
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The  second  method,  which  is  used  in  this  study,  consists  of  evalu- 
ating the  pressure  gradients  normal  to  the  surface  by  applying  the 
^-momentum  equation  at  the  surface.  The  equation  for  p^  is 
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This  equation  is  numerically  modeled  by  using  a one-sided,  second-order 
accurate  finite-difference  algorithm  to  calculate  the  surface  points  of 
interest.  The  surface  pressure  is  calculated  by 


Pn  - — pn  _ i pn  + Q(u,v,w)  AS 
i,j  3 i,j±l  3 i,j±2 


with  1*1.  and  l<i<i  . . The  numerical  representation,  Q,  of 

J Jwing  wing 


v 


the  stress  terms  in  Eq  75  uses  second-order  central  difference  terms 
for  those  graidents  not  on  the  surface. 


For  the  total  delta  wing,  the  same  wing  grid  points  model  both 
the  upper  and  lower  transformed  wing  surfaces.  When  the  MacCormack 
predictor  and  corrector  terms  are  applied  to  the  computational  plane 
below  the  wing,  the  wing  nodal  points  model  the  lower  surface.  When 
the  integration  is  occurring  above  the  wing,  the  wing  modal  points 
model  the  upper  surface  values.  The  leading  edge  is  treated  as  a 
singularity  point  and  is  discussed  in  a later  subsection. 

For  the  cone,  a simplier  boundary  condition  is  used  to  calculate 
the  surface  pressure.  From  boundary  layer  theory,  we  know  that 
Q 2 0,  thus  the  pressure  gradient  normal  to  the  surface  is  zero 
and  the  pressure  pan  be  calculated  by  the  finite-difference  equation 


pn  = A pn  Ipn 

ri,j  3 i,j+l  3 i , j+2 


For  both  bodies,  the  density  at  the  surface  is  determined  by  apply- 
ing the  surface  pressure  and  temperature  in  the  equation  of  state. 


Symmetry  Boundary  Conditions.  A symmetry  boundary  condition 
exists  in  the  transformed  plane  for  both  the  delta  wing  and  the  cone. 


where 
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The  "reflection"  method  is  used  to  calculate  the  primitive  variables  on 
the  plane  of  symmetry.  The  finite-difference  equations  used  for  this 
reflection  technique  are 


i-l.J 


i+l.J 


(79) 


v 


= 0 


and 


*i-l.J  ' *i+l,J 


(80) 


where  i = 2 and  <J>  = u,  w,  p,  p,  and  e.  The  i - 1 column  represents  an 
additional  column  of  grid  points  placed  next  to  the  plane  of  symmetry. 
This  column  of  grid  points  causes  the  integration  along  the  line  of 
symmetry  to  perceive  that  a mirror  image  solution  is  evolving  from  the 
opposite  side  of  the  symmetry  plane.  This  technique  provides  an  accurate 
solution  for  the  flow  field  qualities  on  the  plane  of  symmetry  with  a 
minimal  amount  of  computation. 

Free  Stream  Boundary  Conditions.  The  free  stream  boundary  con- 
ditions are  applied  on  the  extremities  of  the  computational  domain. 

These  conditions  are  applied  under  the  assumption  that  the  free  stream 
or  infinity  boundaries  completely  encompass  the  disturbed  flow  region 
of  the  body.  The  flow  through  these  boundaries  must  be  supersonic  with 
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respect  to  the  transformed  coordinates.  When  this  condition  is 
satisfied,  the  variables  at  the  exterior  grid  points  can  be  set  equal 
to  the  free  stream  conditions  and  held  fixed  during  the  entire  in- 
tegration process.  These  free  stream  conditions  for  uniform  flow  at 
an  angle  of  attack  are 
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(83) 


(84) 


(85) 


e *=  2^  ! (u2  + v2  + w2)  <86) 

where  a is  the  angle  of  attack,  M is  the  free  stream  Mach  number, 

and  V is  the  characteristic  reference  velocity, 
max 


Leading  Edge  Singularity.  The  leading  edge  singularity  point  re- 
quires special  attention.  The  boundary  conditions  at  this  point  are 
identical  to  the  surface  boundary  conditions.  However,  numerical  con- 
siderations require  (Ref  119)  that  the  properties  at  this  singularity 
point  be  triple  valued  in  the  numerical  computation. 


58 


The  no-a lip  conditions  and  the  Isothermal  conditions  arc  applied 

at  the  leading  edge.  The  pressure  and  density  are  triple  valued  at 

the  point.  The  three  values  of  pressure  are  denoted  as  upper,  side, 

and  lower  values  (V  , 1'  , and  1’,)  and  are  assumed  to  exist  slmultane- 
u s I. 

ously  at  the  leading  edge.  These  values  are  calculated  from  the  f.  and 
H momentum  equations.  The  three  values  of  density  are  determined  from 
the  equation  of  state. 

Initial _Condl t Ions 

To  Initiate  the  numerical  Integration,  the  flow  field  quantities 

at  each  grid  point  In  the  computational  plan**  must  he  specified.  It 

has  been  numerically  shown  that  two-step  finite-difference  schemes  may 

be  unstable  If  these  Initial  conditions  differ  greatly  from  the  steady 

state  solution  (Refs  120  and  121).  Since  the  stability  of  linear 

systems  Is  not  affected  by  the  choice  of  Initial  conditions,  then  this 

phenomenon  represents  a nonlinear  Instability  in  the  governing  equations. 

As  a result  of  this  nonlinear  Instability,  the  MacCormack  finite- 

difference  scheme  is  not  absolutely  convergent.  This  means  that  for  an 

(o) 

arbitrary  Initial  vector  X , where  X Is  an  n-eomponent  vector,  the 
(k) 

vector  sequence  X may  not  always  converge.  In  tact,  theorems  which 
prescribe  the  sufficient  conditions  for  convergence  of  nonlinear 
Iterative  methods  (Kef  122)  generally  require  that  the  Initial  vector 
X<‘°  closely  approximate  the  presumed  final  solution.  A,  ot  the  con- 
verged flow  field.  This  "nearness"  Is  usual Iv  expressed  In  terms  of 
the  vector  norm  | ) X ^ ^ - A | | . 

In  tills  Investigation,  the  Initial  starting  conditions  are  speel- 
t t . I it  ,-ach  Interior  grid  point  on  the  comput  at  Iona  1 plane.  These  grid 


points  arc  assigned  values  of  the  free  stream  conditions  given  by 
81-8h.  The  Integration  t lino  stop  Is  determined  from  the  von  Neumann 
stability  orttorla  as  discussed  In  Appoiultx  0. 

Tho  numerical  Integration  Is  Initiated  hv  an  "impnlsivo"  start  In 
whloh  tho  body  ts  Instant anoonsly  accelerated  from  rost  to  t roo  stream 
conditions.  Tho  no-sllp  boundary  conditions  on  t hi'  hodv  surfaoo 
croatos  olthor  a compress  Ion  or  expansion  wave  that  propagates  Into  tho 
flow  field  and  eventual  Iv  sett  It'S  down  as  the  solution  converges.  II 
the  compression  or  expansion  wave  Is  strong,  then  numerical  oscillations 
are  created  which  may  lead  t o pro  (tram  failure.  However,  hv  apply  Inc 
numerical  damping  through  the  addition  o!  an  artlllclal  vlscosltv,  these 
nonltne.'ir  lustahl  Miles  can  he  sutflclentlv  dampened  so  that  a steadv 
state  solution  can  he  achieved. 

An  alternative  method  was  used  when  a series  t'l  ancles  t'l  attack 
for  the  same  I low  conditions  were  investigated.  The  previous  con- 
verged solvit  ton  was  used  as  the  Initial  condition  lor  the  next  angle 
of  attack  case.  Convergence  time  was  reduced  hv  lOt  using  this  method. 

St  ah l 1 1 1 y Ana  I vs  I s 

The  numerical  solution  ot  the  governing  equations  approximates  the 
exact  solution,  II  the  finite  difference  equations  and  tlioli  solid  Ions 
satisfy  the  requirements  ot  stability,  consistency,  and  convergence. 

In  tills  Investigation,  stability  l s directly  associated  with  the  size 
of  the  time  step  In  the  matching  dtiection.  Tho  maximum  time  Increment 
by  which  the  numerical  solution  mav  hi'  advanced  at  any  time  step  Is 
dependent  on  the  spatial  coordinate  gild  sl/.e,  the  location  ol  the 
nodal  grid  points,  and  the  solut  Ion  itself.  tu  order  to  examine  the 
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stability  requirements  of  this  finite-difference  technique,  a definition 
for  stability  Is  needed.  Let  U(n,£,t)  be  the  exact  solution  of  the 
governing  equations  and  let  U(iAn,jA£,nAt)  represent  the  fi  Ite-dlffer- 
ence  solution  of  the  same  equations.  The  error,  E,  of  the  finite-differ- 
ence approximation  is  defined  as  |u(iAr|,  jA£,nAt)  - ll(n,£,t)|.  The 
finite-difference  equations  are  stable  If  the  following  two  conditions 
are  satisfied:  (1)  E remains  bounded  as  n for  a fixed  Ail,  A£,  and  At 
and  (2)  E remains  bounded  as  the  mesh  is  refined  (i.e.,  as  At,  An,  and 
A£-*0)  for  a fixed  value  of  nAt.  Thus,  any  finite-difference  scheme  which 
allows  the  growth  of  the  numerical  error  E with  time,  eventually  "de- 
stroys" the  true  solution  of  the  equations  and  thus  becomes  unstable. 

In  this  study,  the  maximum  time  step  for  numerical  stability  Is 
determined  by  the  Courant-Friedrichs-Lewy  (CFL) , the  diffusion,  and  the 
mixed-derivative  stability  requirements.  By  applying  a von  Neumann 
stability  analysis  separately  to  the  linearized  inviscid,  diffusion, 
and  mixed  derivative  parts  of  the  governing  equations,  a maximum  time 
step  is  developed  (Ref  123).  This  time  step  is 
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where  c is  an  adjustable  constant  and  At,,,,,.  At,,,,,  and  At,.,,,,  are  defined 

INV  l)r  MXD 

in  Appendix  C.  A complete  development  of  this  stability  criteria  is 
given  in  Appendices  C and  D. 

Several  investigators  (Refs  122-123),  using  the  shock  capturing 
method,  have  encountered  numerical  oscillations  in  the  vlclnltv  of 
strong  shock  waves.  These  nonlinear  Instabilities  can  be  overcome  by 
adding  an  artificial  viscosity.  In  this  study,  two  types  of  artificial 
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viscosity  terms  are  used  to  dampen  oscillations.  A normal  stress  damp- 
ing terra  by  McRae  (Ref  94)  is  used  for  initial  transients  and  a fourth- 
order  damping  term  by  MacCormack  Is  used  for  shocks. 

During  program  startup,  large  impulsive  oscillations  are  developed 
near  the  body  surface.  These  initial  transients  can  be  large  enough  to 
cause  numerical  instability.  By  applying  a multiplier  to  the  second 
coefficient  of  viscosity,  these  oscillations  are  effectively  reduced. 

As  the  solution  progresses  in  time,  this  multiplier  is  gradually  re- 
duced to  its  minimum  value  of  one.  This  numerical  damping  technique 
adversely  affects  the  maximum  time  step  that  can  be  used,  as  seen  in 
Appendix  C. 

For  numerical  oscillations  in  the  vicinity  of  strong  pressure 
gradients,  a fourth-order  damping  term,  suggested  by  Tannehill  et  al 
(Ref  133),  is  adopted  with  modifications.  Instead  of  using  a rigorous 
transformation  of  the  second-order  derivatives  of  density  into  the 
transformed  space,  the  second-order  derivatives  are  approximated  by 
the  derivatives  with  respect  to  the  transformed  Independent  variables 
( r| » C)  - The  net  results  are  two  artificial  viscosity-like  terms  of  the 
form 

Di,j  “ X C«.  ( K£+lV’  ~ *W'  ) (88) 
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where  the  U's  are  defined  by  Kq  30,  the  summations  on  £ indicate  one 
term  for  each  of  the  two  spatial  directions  (i.J),  the  c's  are  constant 
coefficients,  the  K's  are  variable  coefficients  defined  by 
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and  the  6^(  ) and  A^(  ) operators  are  forward  and  backward  dif- 

ferences, respectively,  defined  by 
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When  smoothing  is  desired  D,  . is  added  to  Eq  57  on  the  predictor  step 

i > J 

and  D?^T  is  added  inside  the  outer  brackets  of  Eq  58  on  the  corrector 

i t J 


The  variable  coefficients  are  composed  of  a normalized  second- 

order  difference  of  flow  field  density.  The  subscript  indicates 
which  direction  the  difference  is  in  (i  or  j)  and  also  the  center  of 
the  difference.  The  superscript  Indicates  the  time  level  used  for 
the  values  of  density.  This  variable  coefficient  is  essentially  zero 
in  smooth  regions  of  the  flow  field  and  approaches  a maximum  value  of 
one  in  regions  of  large  point-wise  oscillations.  The  theoretical 
maximum  value  of  the  coefficient  product  c^K^  is  less  than  or  equal  to 
0.5  to  prevent  this  term  from  causing  an  instability  itself.  The 
coefficients  can  theoretically  reach  a value  of  one,  causing  c^  co- 
efficients to  be  restricted  to  a value  of  one-half  or  less.  However, 
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since  coefficients  are  usually  much  smaller  than  one,  the  c^  co- 
efficients can  be  larger  than  one-half.  Smaller  amounts  of  damping  or 
damping  in  only  one  of  the  two  spatial  directions  can  also  be  used. 

This  is  achieved  by  independently  setting  the  constant  coefficients  in 
each  of  the  two  directions  equal  to  the  appropriate  values. 

Consistency  and  Order  of  Accuracy 

In  order  ta  examine  the  consistency  and  order  of  accuracy  of  a 
finite-difference  scheme,  it  is  necessary  to  understand  what  partial 
differential  equation  is,  in  fact,  being  solved  by  the  finite-difference 
algorithm.  This  differential  equation  Is  called  the  modified  equation 
and  aside  from  the  round-off  error,  actually  represents  the  original 
partial  differential  equation  when  a numerical  solution  is  computed. 

The  modified  equation  is  derived  by  first  expanding  each  term  of  the 
finite-difference  equation  in  a Taylor  series  and  then  eliminating 
time  derivatives  of  higher  than  first-order  (including  mixed  time  and 
space  derivatives)  by  the  algebraic  manipulations  described  by  Warming 
and  Hyett  (Ref  124).  The  terms  appearing  in  the  modified  equation  which 
are  not  in  the  original  partial  differential  equation  represent  a form 
of  truncation  error  introduced  by  the  finite-difference  scheme.  These 
error  terms  allow  one  to  determine  the  order  of  accuracy  and  consistency 
of  the  finite-difference  approximation. 

If  we  examine  the  linearized  form  of  the  governing  equation,  we 
find  that  the  modified  equation  is 
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where  the  matrices  A,B,C,D,E,  and  F are  defined  in  Appendix  C.  The 
order  of  accuracy  of  the  finite-difference  scheme  is  defined  by  the 
lowest-order  powers  of  the  increments  At,  An,  and  AC  appearing  in  the 
error  terras.  Thus,  according  to  the  modified  equation,  the  MacCormack 
scheme  is  a second-order  accurate  algorithm  for  the  linear  governing 
equations.  This  definition  of  order  of  accuracy  is  consistent  with  that 
used  by  Richtmyer  and  Morton  (Ref  109)  in  which  truncation  error  is 
normally  examined. 

A finite-difference  scheme  is  defined  as  consistent  if  the  dif- 
ference between  the  partial  differential  equation  (with  the  initial 
and  boundary  conditions)  and  the  finite-difference  equation  tends  to 
zero  as  At,  Ar|,  and  AC  approach  zero  in  some  arbitrarily  stable  manner. 
This  implies  that  the  error  terras  in  the  modified  equation  approach 
zero  in  the  limit  as  At,  An,  and  AC"*0.  By  examining  Eq  94,  it  is 
obvious  that  for  the  MacCormack  scheme,  the  finite-difference  equation 
is  consistent  with  the  original  partial  differential  equation.  I f by 
some  special  relationship  the  error  terms  did  not  tend  to  zero  as  the 
increments  vanish,  then  the  finite-difference  scheme  is  said  to  be  in- 
consistent. In  this  situation,  the  modified  equation  yields  directly 
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the  differential  equation  with  which  the  finite-difference  scheme  is 
consistent. 


Convergence 

In  order  to  show  equivalance  between  the  finite-difference  equations 
and  the  partial  differential  equations,  it  is  necessary  to  demonstrate 
that  the  MacCormack  scheme  is  convergent.  Richtmyer  and  Morton  (Ref 
109)  define  a convergent  finite-difference  scheme  as  one  in  which  the 
numerical  solution  approaches  the  solution  of  the  partial  differential 
equation  as  the  step  size  approaches  zero.  The  proof  of  convergence  is 
particularly  complex  for  the  finite-difference  scheme  associated  with 
the  governing  equations.  This  is  especially  true  when  boundary  con- 
ditions are  included  in  the  analysis.  Thus,  in  this  investigation,  only 
the  linearized  form  of  the  governing  equations  is  examined. 

According  to  Lax's  Equivalence  Theorem  for  an  initial  value  problem, 
a finite-difference  solution  converges  to  the  exact  solution  if  and  only 
if  1)  the  differential  problem  is  properly  posed  and  therefore  possesses 
a genuine  solution,  2)  the  finite-difference  equations  are  consistent, 
and  3)  the  numerical  computation  is  stable.  In  this  analysis,  it  is 
assumed  that  the  initial  value  problem  is  properly  posed.  Tills  implies 
that  a solution  of  the  governing  equations  exists,  is  unique,  and  depends 
continuously  on  the  initial  data.  Therefore,  by  satisfying  the  criteria 
of  Lax's  Equivalence  Theorem,  it  is  assumed  that  the  MacCormack  second- 
order  finite-difference  algorithm  is  convergent  for  the  governing 
equations . 

The  numerical  analysis  of  the  MacCormack  method  has  been  based  on 
the  mathematical  theory  of  linear  partial  differential  equations.  By 


examining  the  linearized  form  of  the  governing  equations,  a heuristic 
assessment  has  been  made  of  the  MacCormack  algorithm  as  applied  to  the 
governing  equations.  Although  there  is  no  guarantee  that  the  linear 
and  nonlinear  characteristics  of  the  finite-difference  scheme  are 
similar,  there  is  ample  evidence  which  tends  to  support  this  assumption. 
Thus,  in  the  absence  of  a systematic  analysis  the  numerical  schemes 
associated  with  the  governing  equations,  it  is  assumed  that  the  above 
analysis  will  provide  a practical  insight  into  the  characteristics  of 
the  MacCormack  finite-difference  scheme. 
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IV.  Numerical  Computation 

In  order  to  obtain  a numerical  solution  for  the  supersonic  and 
hypersonic  flow  around  a delta  wing,  it  is  necessary  to  translate  the 
finite-difference  equations  and  the  initial  and  boundary  conditions 
into  a computer  code.  The  purpose  of  this  chapter  is  to  describe,  in 
general  terms,  the  basic  structure  and  operation  of  the  computer  codes 
used  in  this  investigation.  These  computer  codes  use  the  MacCormack 
finite-difference  scheme  to  solve  the  conical,  viscous  governing 
equations.  The  input  data,  sequence  of  operations,  and  the  calculated 
results  of  these  codes  are  discussed  in  this  chapter.  A brief  descrip- 
tion of  each  subroutine  used  in  the  computer  programs  is  given  in 
Appendix  E . 

Computer  Codes 

Three  basic  computer  codes  are  used  in  this  study.  These  codes 
are  the  CONE,  DELTA1,  and  DELTA  computer  programs.  They  are  used  to 
solve  the  flows  over  a circular  cone,  flat  delta  wing,  and  around  a 
thin  planar  delta  wing,  respectively.  The  basic  structure  and  opera- 
tion of  each  of  these  programs  are  the  same.  The  differences  in  these 
codes  occur  in  how  the  boundary  conditions  are  modeled  and  in  what 
coordinate  transformations  are  used.  Because  of  the  close  similarity 
among  the  programs,  only  one  code  is  described  in  this  chapter.  This 
is  the  DELTA  computer  program.  Differences  between  the  basic  computer 
programs  used  in  this  study  are  highlighted  in  the  description  of  the 
DELTA  program. 

DELTA 

The  computer  code,  DELTA,  was  written  to  solve  the  generalized 
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conical  governing  equations  numerically  and  to  validate  the  conical 
approximation  technique.  With  only  minor  changes  in  the  input  data  and 
internal  code  structure,  the  code  is  designed  to  solve  any  problem  which 
satisfies  the  conical  flow  approximation.  The  program  Is  structured  in 
a modular  or  subprogram  form  so  that  each  module  or  subprogram  can  be 
coded  and  debugged  separately.  This  procedure  adds  to  the  total  com- 
putation time  of  the  program  but  it  minimizes  the  number  of  changes 
needed  to  solve  various  conical  flow  problems  and  it  reduces  the  time 
required  to  debug  the  program.  No  major  effort  was  made  to  reduce  the 
running  time  or  storage  requirements  of  the  computer  program.  The 
primary  emphasis,  in  developing  this  code,  was  to  keep  the  program  as 
flexible  and  generalized  as  possible  so  that  various  changes  in  boundary 
conditions  and  coordinate  transformations  could  be  easily  made.  This 
code  represents  the  first  step  towards  development  of  a user-oriented 
code  which  can  be  used  to  solve  a variety  of  conical,  viscous  flow  field 
problems . 

The  DF.LTA  program  begins  its  calculation  by  initializing  the  grid 
points  in  the  computational  plane  with  free  stream  and  boundary  con- 
dition values.  A constant  time  step  is  specified  at  the  beginning  of 
each  run  and  a maximum  number  of  iterations  is  given  as  part  of  the 
input  data.  The  computer  program  is  run  until  it  reaches  the  maximum 
number  of  iterations  or  satisfies  the  convergence  criteria.  When 
either  of  these  two  conditions  are  reached,  the  calculations  are 
terminated  and  the  results  are  printed  out  and  stored  on  magnetic  tape. 
If  an  interim  solution  is  stored  on  magnetic  tape,  then  this  solution 
is  used  as  the  initializing  data  for  the  next  computer  run  and  the 
calculations  are  continued.  When  the  convergence  criteria  is  finally 


I 


69 


satisfied,  the  numerical  integration  is  terminated  and  the  final  results 


are  printed  out  and  stored  on  magnetic  tape.  These  final  results  in- 
clude: 

1.  The  nondimens ional  scalar  velocity  components  in  the  x,  y,  and 
z directions. 

2.  The  normalized  pressure,  density  and  internal  energy  as  defined 
in  chapter  2. 

3.  The  cross-flow  Mach  number,  M , and  the  slope  of  the  conical 
streamlines,  Yc»  with  respect  to  the  n axis  as  given  by 
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These  quantities  are  calculated  for  each  point  in  the  computational 
plane. 

No  capability  for  interpolating  between  the  computational  grid 
nodes  is  provided  in  the  program.  Therefore,  it  is  necessary  to  insure 
that  a spatial  grid  point  is  specified  at  each  location  for  which  flow 
field  data  is  required. 


Sequence  of  Calculations 

In  order  to  initiate  the  numerical  integration,  the  DELTA  program 
is  required  to  receive: 

1.  The  free  stream  Mach  number  and  angle  of  attack.  The  angle  of 
attack  is  measured  relative  to  the  body  coordinate  system  as  shown  in 
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Figures  8a.  and  9a. 

2.  The  sweep  angle  for  the  DELTA  and  DELTA1  programs  and  the  cone 
half  angle  for  the  CONE  program. 

3.  The  free  stream  isentropic  total  temperature  and  the  isothermal 
surface  temperature,  in  degrees  Rankine. 

4.  The  free  stream  Reynolds  number  and  a reference  length.  This 
reference  length  is  defined  as  the  distance  from  the  vertex  to  the 
desired  y-z  plane. 

5.  The  time  step  size  and  the  maximum  number  of  time  step  itera- 
tions. 

6.  The  number  of  grid  points  on  the  body  surface  and  in  the  com- 
putational plane.  The  location  of  these  grid  points  is  determined  by 
inputing  A£  and  by  calculating  An. 

7.  The  appropriate  coordinate  transformation,  as  defined  in 
Appendix  A,  and  the  grid  location  of  the  n-axis  for  the  DELTA  program. 

8.  The  normal  stress  damping  coefficients  and  the  pressure  damping 
coefficients. 

Once  the  above  input  data  is  supplied,  the  DELTA  program  begins  its 
calculations  in  the  following  order: 

1.  The  computational  grid  points  are  initialized  with  free  stream 
and  boundary  condition  values  or  with  an  interim  flow  solution,  which- 
ever is  appropriate. 

2.  The  coordinate  transformation  derivatives  are  calculated  next 
in  the  COORD  subroutine. 

3.  At  this  point,  the  main  program  DO  loop  is  entered  for  the 
repetitive  calculation  of  all  the  flow  field  quantities  at  each  time 
step  until  convergence  or  the  maximum  number  of  time  step  iterations  is 
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reached . 


4.  Within  this  DO  loop,  the  local  Reynolds  number  is  calculated  for 
each  grid  point. 

5.  The  predictor  term  of  the  MacCormack  finite-difference  scheme 
is  calculated  next  by  calling  the  subroutine  PREDICT.  In  this  sub- 
routine, the  surface  boundary  conditions  are  computed  first  by  calling 
subroutine  BOUND.  After  that,  a double  DO  loop  is  entered  wherein  the 


following  calculations  are  made  for  each  grid  point  (i,j) : 

(a)  The  pressure  damping  terms  (D  n.)  and  (D.n^)  are  deter- 

i » J 1 > J 

mined  by  calling  the  subroutines  DAMPF  and  DAMPG. 

(b)  The  subroutine  VECTOR  is  called  next  to  solve  for  the 


stress  and  heat  conducting  terms  as  well  as  determine  the  values  of  the 
F,  G,  and  H vectors. 

(c)  The  predicted  value  of  the  j vector  at  the  new  time 
level  is  computed  by  calling  subroutine  DECODE. 

(d)  The  subroutine  SOLVE  is  used  next  to  solve  for  the  flow 
quantities  from  the  calculated  ^ vector  at  the  new  time  level. 

6.  The  flow  qu...tities  on  the  opposite  side  of  the  symmetry  plane 
are  determined  by  calling  subroutine  SYM. 

7.  The  local  Reynolds  number  is  recalculated  with  the  flow 

quantities  from  the  U,  . vector  at  the  new  time  level. 

i » J 

8.  The  corrector  term  of  the  MacCormack  finite-difference  scheme 
is  calculated  next  by  calling  subroutine  CORRECT.  This  subroutine 
duplicates  the  predictor  subroutine  (5)  except  for  some  changes  in 
the  finite-difference  representatives  and  flow  quantity  levels.  The 
same  subroutines  and  sequence  of  events  are  followed. 


9.  After  the  corrector  step,  the  SYM  subroutine  is  called  to 
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recalculate  the  flow  quantities  on  the  opposite  side  of  the  symmetry 
plane . 

10.  The  final  step  in  the  main  DO  loop  is  to  test  the  new  solu- 
tion for  convergence.  If  the  new  solution  is  not  converged,  the 
numerical  integration  process  is  begun  again  at  step  (4)  until  the  con- 
vergence criteria  is  met  or  the  maximum  number  of  time  steps  is  reached. 

11.  At  the  end  of  the  numerical  integration  cycle,  the  flow  field 
results  are  printed  out  and  the  interim  or  final  solution  is  stored  on 
magnetic  tape. 

The  computational  sequence  of  this  code  is  shown  in  Figures  11  and 

12.  This  program  is  repetitively  run  until  a converged  solution  is 


obtained. 


CALL  BOUND 


Fig.  12.  Schematic  Diagram  of  PREDICT  and  CORRECT 
Subroutines. 
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V.  Numerical  Results 


In  this  chapter,  the  numerical  solutions  for  both  the  supersonic 
flow  around  a cone  and  the  supersonic/hypersonic  flow  around  a thin 
planar  delta  wing  are  discussed.  In  each  of  these  cases,  the  conical, 
viscous  approximation  is  applied  to  the  governing  equations  to  determine 
the  laminar  flow  field  characteristics.  The  numerical  results  of  these 
calculations  are  compared  with  various  analytical  solutions  and  ex- 
perimental data. 


Cone  Flow  Analysis 

In  order  to  verify  the  accuracy  of  the  computer  codes,  the  super- 
sonic flow  around  a cone  at  zero  angle  of  attack  was  computed  and  com- 
pared with  previous  analytical  and  experimental  results.  The  coordinate 
transformation  used  in  this  calculation  is  discussed  in  Appendix  A and 
the  physical  and  transformed  planes  are  shown  in  Figures  8a  and  8b.  The 
free  stream  and  surface  boundary  conditions  chosen  for  this  calculation 
are 
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where  the  Reynolds  number  is  based  on  a reference  length  of  x=4.0  inches. 
These  flow  conditions  arc  identical  to  those  used  by  Tracy  (Ref  127)  in 
his  experimental  studies  of  hypersonic  flow  over  a yawed  circular  cone. 
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In  this  investigation,  two  different  rectangular  grid  systems  were 
used  to  calculate  the  flow  field.  These  included  a coarse  mesh  with  a 
radial  step  size  of  A£»0.0079  and  a fine  mesh  with  a step  size  of 
A£>*0. 00395.  A circumferential  grid  size  of  Ari"3°  was  used  in  both 
systems.  The  computational  domains  contained  63  equally  spaced  n grid 
points  ranging  from  -3°<r)£l830 , and  30  constantly  spaced  £ grid  points 
from  0^-10°  to  the  free  stream  boundary  surface.  The  free  stream 
boundaries  were  located  at  £-tan(16.21°)  (fine  mesh)  and  at  £-tan(22.06°) 
(coarse  mesh) . Both  of  these  boundary  locations  were  far  enough  from 
the  cone  surface  and  bow  shock  so  as  not  to  affect  the  numerical  solu- 
tion. 

Figures  13-16  illustrate  a comparison  of  the  numerical  results  of 
this  technique  with  the  solutions  obtained  by  McRae.  The  squares  and 
triangles,  in  these  figures,  represent  the  numerical  results  for 
A£»0. 00395  and  AC”0.0079,  respectively.  The  solid  curve  depicts  the 
numerical  solution  obtained  by  McRae.  In  all  of  these  solutions,  no 
numerical  damping  was  used.  The  undamped  solutions  were  obtained  in 
order  to  compare  the  results  more  easily  and  to  assess  the  accuracy  of 
the  computer  code. 

The  static  pressure  distribution  in  the  0 direction  above  the  cone 

surface  is  shown  in  Figure  13.  The  normal  pressure  gradient  in  the 

boundary  layer  is  zero  and  the  numerical  solutions  near  the  cone  surface 

are  within  one  percent  of  each  other.  If  we  assume  that  the  attached 

bow  shock  location  coincides  with  the  static  pressure  transition 

centerpoint  (a  definition  used  by  Tracy),  then  the  experimental  shock 

location  is  at  0-0  "3.56°  (Ref  127).  The  numerical  shock  locations  for 
c 

A£»0. 00  395  and  AF.-0.0079  are  at  O-O  *3.57°  and  O-O  -3.52°,  respectively. 
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Fig.  13.  Comparison  of  Static  Pressure  vs.  0 Above  the  l-one  Surtaee 
for  Otfferent  CrUl  Sizes  with  McRae's  Solution. 
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McRae's  solution  shows  the  shook  at  3.48°  above  the  cone  surface.  These 
results  Indicate  that  the  flue  mesh  solution  agrees  slightly  better  with 
experimental  data  than  does  the  coarser  mesh  solution.  However,  both 
numerical  results  yield  excellent  agreement  with  both  McRae's  solution 
and  Tracy's  experimental  data. 

The  surface  pressure  on  the  cone  la  determined  from  Kq  77,  where 
the  normal  pressure  gradient  Is  zero . The  nond l mens  Iona  1 Iced  surface 

* >4 

pressures  for  both  the  line  and  coarse  grid  systems  are  4. ‘>0x10  and 
4.54x10  , respectively.  The  surface  pressure  calculated  bv  McRae  Is 
4.53x10  . Experimental  results  by  Tracy  show  the  normalized  surface 
pressure  to  he  4.t>JxlO  ' . This  discrepancy  between  calculated  and 
measured  results  can  be  attributed  to  experimental  error.  Instrumenta- 
tion error,  and  numerical  modeling  error.  It  can  be  seen  that  the 
simplified  model  of  the  surface  boundary  conditions  provides  excellent 
agreement  with  McRae's  solution  and  good  agreement  with  Tracy's  ex- 
perimental data. 

In  Figure  14,  the  total  velocity  distribution  In  the  0 direction 
is  presented  for  x-*-».0  Inches.  The  numerical  resolution  of  the  boundary 
layer  and  bow  shock  Is  marginal  tor  the  tine  grid  system  and  is  poor  for 
the  coarser  grid  system.  However,  the  agreement  between  McRae's  cal- 
culations and  the  numerical  results  Is  remarkable.  In  general,  there 
is  little  difference  between  the  fine  and  coarse  mesh  systems,  except 
that  the  fine  mesh  solut Ion  provides  a much  better  definition  ot  the 
flow  field.  Both  numerical  results  yield  nearly  identical  behavior  as 
that  noted  In  similar  experimental  observat Ions  bv  Tracv. 

The  static  temperature  distribution  In  the  0 direction  above  the 

3T 

cone  surface  Is  shown  In  Figure  15.  It  can  be  seen  that  0 and 
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cone  surface.  Tills  Is  occurring  because  of  the  large  amount  of  friction- 
al heat  being  generated  In  the  boundary  layer,  even  though  T ''T^. 

Because  of  the  limited  number  of  grid  points  in  the  viscous  region,  a 
verification  of  this  result  was  made  with  a .-.ero-pressure  gradient 
similarity  solution  of  the  boundary  layer  equation.  In  order  to  convert 
the  similar  solution  results  given  by  l.ow  (Ref  libl  from  the  non- 
dimensional  Blasius  variables  to  the  present  nondtmenslonal  variables. 

It  was  necessary  to  specify  V,  p , and  T at  the  boundary  layer  edge  as 
well  as  at  the  cone  surface.  Rather  than  choose  free  stream  conditions 
for  the  outer  edge  conditions  as  would  be  consistent  with  classical 
boundary  layer  theory,  the  viscous  edge  conditions  were  chosen  from  the 
fine  mesh  solution.  A survey  ot  calculated  impact  pressures  in  the  0 
direction,  as  showft  In  Vigure  In,  was  used  to  determine  the  boundary 

laver  edge  (0  -0  The  normalised  surface  pressure  p « 4. SO  x 

* e c 

-4 

10  was  assumed  tv'  be  constant  throughout  the  viscous  region.  This  pres- 
sure coupled  with  the  edge  temperature  T /T  - 0.117  gave  the  edge 

' »'  O 

•> 

density  p /p  - t.84  x 10  '.  From  Figures  14  and  lb,  the  edge  velocity 

e o 

was  found  tv’  ho  V /V  ■ 0.4.,.  Rv  applying  these  edge  condit  ions  tv’  the 
e max 

boundary  laver  equations,  a similar  temperature  profile  was  obtained, 
thus  confirming  the  correctness  of  the  numerical  results. 

The  effect  of  grid  refinement  in  the  £ direction  is  most  signifi- 
cant in  the  vicinity  ot  the  shock  wave,  as  shown  In  Figures  Id  and  If. 

For  As.“0.00»5  and  A0-O.22f°  (McRae’s  step  sire),  the  numerical  oscilla- 
tions are  moderate  on  both  sides  of  the  shock.  However,  when  .V.  0.0074, 
the  oscillations  on  the  upwind  si  vie  of  the  shock  become  larger.  This 
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causes  the  numerical  instabilities  (Gibbs  phenomenon)  to  be  spread  over 
a larger  distance  before  they  are  completely  damped  out. 

In  conclusion.  It  can  be  seen  that  this  computer  code  accurately 
calculates  the  supersonic,  viscous  flow  field  around  a cone.  The 
velocity,  temperature,  and  pressure  results  compare  quite  favorably  with 
both  McRae's  solution  and  Tracy's  experimental  data.  The  fine  mesh 
solution  provides  just  enough  flow  field  definition  in  order  to  deter- 
mine the  flow  characteristics  adequately  in  the  boundary  layer  and  near 
the  bow  shock. 

Delta  W ins  Expans  ion  S i do  FI ow  Analysis 

The  next  phase  in  this  investigation  is  to  examine  both  the  super- 
sonic and  hypersonic  flows  over  the  leeside  of  a planar  delta  wing.  The 
purpose  of  this  analysis  is  to  verify  the  applicability  of  the  conical, 
viscous  approximation  in  solving  the  expansion  side  flow  field.  In 

) 

order  to  do  this,  two  simplifications  were  made  in  modeling  the  flow. 
These  were  (1)  to  assume  the  flows  above  and  below  the  wing  surface  are 
independent  of  each  other  and  (2)  to  neglect  some  of  the  details  of  the 
complex  flow  in  the  immediate  vicinity  of  the  delta  wing  vertex  and 
leading  edge. 

In  the  first  assumption,  the  upper  and  lower  wing  flows  are  in- 
dependent of  each  other  because  the  leading  edge  flow  is  supersonic. 

This  follows  from  the  fact  that  in  supersonic  linear  theory,  disturbances 
are  propagated  downstream  in  a Mach  cone  which  has  its  axis  aligned  with 
the  free  stream  and  which  has  a semi-angle  equal  to  the  free  stream  Mach 
angle.  Thus,  any  disturbance  generated  by  either  wing  surface  is  con- 
fined to  a Mach  cone  which  does  not  intersect  the  leading  edge,  since 
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the  edge  is  less  swept  than  the  Mach  cone.  However,  in  a real  flow, 
this  does  not  hold  since  the  regions  of  influence  are  determined  by 
local,  rather  than  free  stream  conditions  (Ref  128).  For  delta  wings, 
the  two  surfaces  are  independent  if  the  bow  shock  is  attached  at  the 
leading  edge.  Wing  thickness  and  viscous  development  effects  can 
detach  the  shock  wave  and  thus  cause  a strong  overflow  from  the  com- 
pression side  of  the  wing.  This  overflow  influences  the  supersonic 
and  hypersonic  flow  on  the  expansion  side  of  a delta  wing,  as  described 
in  Chapter  1.  However,  in  this  numerical  model  the  bow  shock  is  attach- 
ed to  the  leading  edge  and  thus  the  compression  side  flow  does  not  in- 
fluence the  upper  surface  flow  field  characteristics. 

The  second  simplification,  which  is  used  throughout  this  study,  is 
to  neglect  the  three-dimensional  flow  field  effects  near  the  wing  vertex 
and  around  the  leading  edges.  The  flow  over  the  wing  vertex,  at  angle 
of  attack,  is  characterized  by  a bow  shock,  a leading  edge  Prandt 1-Meyer 
expansion,  a pair  of  internal  shock  waves,  and  a developing  boundary 
layer.  These  flow  properties  are  highly  dependent  on  a number  of  param- 
eters, such  as  (1)  leading  edge  geometry,  (2)  angle  of  attack,  (3)  Mach 
number,  (4)  Reynolds  number,  and  (5)  hypersonic  similarity  parameter. 
This  complex  flow  induces  vortex  generation  downstream  of  the  apex 
region,  and  the  development  of  a three-dimensional  boundary  layer. 

Vortex  generation  in  the  boundary  layer  depends  critically  on  the 
cross-flow  and  vlscous-invlscid  interaction  in  the  apex  region,  where 
the  spanwise  pressure  gradients  are  strongest  (Ref  2) . This  has  been 
demonstrated  by  Whitehead  and  Bertram  (Ref  129)  and  by  Rao  (Ref  130). 

The  former  reduced  the  cross- flow  in  the  apex  region  by  rounding  off 
the  sharp  apex,  thus,  locally  "unsweeping"  the  leading  edge.  The 
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latter  drooped  the  portion  of  the  wing  near  the  apex  In  order  to  align 
it  with  the  free  stream  at  angle  of  attack,  thus  making  the  local  angle 
of  attack  in  the  apex  region  zero.  In  both  cases,  the  critical  cross- 
flow  development  In  the  apex  region  was  reduced  which  effectively  sup- 
pressed embedded  vortex  generation.  In  addition,  experimental  measure- 
ments have  shown  that  embedded  vortices  occur  when  \*0(0.1)  (Reis  18 
and  131),  but  not  at  large  values  of  x»  such  as  (Kefs  17  and  25). 

The  reduced  spanwlse  pressure  gradient  (due  to  diminished  expansion  at 
the  leading  edge)  has  the  effect  of  retarding  vortex  generation.  These 
embedded  vortices  are  different  from  those  formed  by  shock-induced 
boundary  layer  separation  at  higher  angles  of  attack  (a>‘>1  ) , as  seen 
In  Cross'  study. 

Rao  and  Whitehead  (Ref  111)  examined  the  hypersonic  boundary  layer 
along  the  centerline  of  a 75°  delta  wing  at  a- 5°  In  M -6.8  flow  (Re/ft  « 
1.2  x 1 0^* ) . It  was  found  that  the  boundary  layer  thickness,  6,  Ini- 
tially develops  as  a two-dimensional  boundary  layer;  however,  a maximum 
d <5 

6 is  reached  and  , becomes  negative.  This  decrease  In  boundary  layer 
dx 

thickness  Is  due  to  the  entrainment  of  low  momentum  fluid  by  the  presence 

of  vortices  in  the  boundary  layer.  Further  downstream,  tills  "trough"  or 

decrease  in  6 begins  to  fill  in  between  the  vortices  as  these  vortices 

move  apart.  An  Increase  In  unit  Reynolds  number  (to  2.0  million  per 

foot  and  to  3.5  million  per  foot)  moves  the  cS  upstream. 

' min 

In  conclusion,  this  investigation  assumes  that  the  effects  ot  the 
viscous  interaction  around  the  leading  edge  and  near  the  wing  vertex  can 
be  neglected  because  those  flows  have  only  a small  effect  on  the  total 
flow  field.  McRae  (Ref  *)/()  and  llankey  (Ref  1 32)  showed  that  the  conical, 
viscous  approximation  provides  good  agreement  with  experimental  results, 
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when  viscous  Interactions  are  weak  (j(<0(l)).  Thus,  the  primary  area  of 
Interest  in  this  research  program  is  to  calculate  the  conical  flow  in 
the  weak-interact ion  region. 


Supersonic  Flow.  In  this  analysis,  only  one  supersonic  case  is  com- 
puted for  the  expansion  side  flow  field.  The  free  stream  conditions 
chosen  for  this  calculation  are 
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where  the  Reynolds  number  is  based  on  a root  chord  length  of  0.173  ft. 
These  flow  conditions  correspond  Identically  to  those  used  by  Bannink 
and  Nebbeling  (Ref  15)  in  their  experimental  Investigation.  As  part 
of  their  study,  they  used  a delta  wing  model  with  a flat  upper  surface 
and  a semi  apex  angle  of  45.3°.  The  wedge  shaped  cross-section  per- 
pendicular to  the  leading  edge  has  an  apex  angle  of  7.5°.  A conical 
head  probe  was  used  to  measure  pitot  pressure  and  flow  angularity. 

The  delta  wing  coordinate  transformation,  as  described  in  Appendix 
A,  was  used  in  these  numerical  calculations  as  well  as  all  other  delta 
wing  calculations.  The  computational  domain,  on  which  these  calcula- 
tions were  made,  consisted  of  a 26(n)x30(O  grid  array,  similar  to  the 
upper  half  of  the  grid  system  shown  in  Figure  9b.  For  the  sweep  angle 
of  45.3°,  the  r)  step  size  was  0.063158  with  14  grid  points  on  the  wing 
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surface.  In  the  £ direction,  a constant,  but  different,  step  size  was 
used  for  each  measurement  height  (pitot  pressure  measurements  by 
Bannink,  Ref  15)  so  that  experimental  and  numerical  results  could  be 
compared  without  interpolation.  These  constant  step  sizes  ranged  from 
0.02393  to  0.02775  for  different  calculations.  The  results  of  these 
calculations  are  shown  in  Figures  17-24. 

In  Figures  17  and  18,  several  spanwlse  pitot  pressure  distributions 
are  shown  for  various  heights  above  the  delta  wing.  The  conical,  viscous 
results  of  tills  technique  are  compared  with  the  inviscld  calculations  by 
Kutler  (Ref  89)  and  the  experimental  measurements  by  Bannink  and  Nebbel- 
ing  (Ref  15).  For  £=0.0718  (Fig  17),  the  experimental  pitot  pressure 
profile  across  the  inboard  shock  is  clearly  affected  by  the  interference 
with  the  wing  boundary  layer.  At  this  height,  the  measured  shock 
strength  is  slightly  less  than  that  measured  at  higher  values  of  £ . This 
numerical  technique  accurately  predicts  the  shock  wave-boundary  layer 
interaction,  which  is  not  accounted  for  in  Kutler's  solution.  Roth 
theoretical  methods  provide  good  agreement  with  experimental  data, 
except  in  the  vicinity  of  the  bow  shock.  In  this  region,  the  current 
numerical  model  does  not  account  for  the  influence  of  the  compression 
side  bow  shock  near  the  leading  edge.  By  neglecting  the  lower  surface 
flow  influence,  a much  weaker  leading  edge  compression  wave  is  computed 
due  to  displacement  effects  of  the  boundary  layer.  This  very  weak  shock 
wave  is  clearly  depicted  in  the  contour  plots  shown  later  in  Figures  22 
and  24. 

Since  the  inviscld  flow  field  and  shock  wave  structure  are  nearly 
conical,  the  velocity  components  in  a spherical  coordinate  system  are 
used  to  delineate  the  three-dimensional  flow  field.  The  velocity  vectors 
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Fig.  18.  Pitot  Pressure  Distribution  in  Spanwise  Direction 
for  Supersonic  Flow  Above  a Planar  Delta  Wing. 
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normal  to  the  rays  through  the  vertex  of  the  delta  wing  are  projected 
on  the  £-n  plane  as  shown  in  Figure  19.  The  relative  magnitude  of 
these  vectors  Is  Illustrated  for  all  but  the  lowest  momentum  region. 
The  direction  of  these  vectors  Is  given  by 


tan  y 

c 


u£-w 

un-v 


(97) 


where  y Is  the  angle  with  the  positive  rj  axis  (see  Appendix  F) . The 
c 

locus  of  these  velocity  vectors  trace  out  "pseudo"  streamlines  which 
converge  at  a vortical  singularity  point  (cross-flow  stagnation  point) 
near  the  origin.  In  the  unperturbed  flow  region  (upper  and  right 
portions  of  Fig  19),  the  conical  cross-flow  velocity  vectors  point  to 
the  origin  of  the  coordinate  system.  Figure  20  illustrates  some  of  the 
cross-flow  conical  streamlines  (780  points).  A linear  interpolation 
technique,  described  in  Appendix  F,  is  used  to  determine  these  steady 
state  path  lines.  Figure  20  shows  more  clearly  how  most  of  these 
conical  streamlines  converge  at  the  vortical  singularity  point,  even 
though  the  numerical  resolution  is  marginal.  This  result  is  consistent 
with  the  experimental  observations  by  Bannink  and  Nebbeling. 

In  Figure  21,  the  calculated  position  of  the  internal  shock  wave 
and  conical  sonic  line  (M  =1.0)  is  shown  in  the  computational  plane. 
Superimposed  on  this  plot  are  the  experimental  results  obtained  by 
Bannink  and  Nebbeling  (Kef  15).  The  measurement  of  the  pressure  rise 
across  a shock  wave  is  dependent  on  the  shock  strength  and  dimensions 
of  the  probe.  Bannink  and  Nebbeling  (Kef  133)  found  that  for  cylindri- 
cal pitot  probes  Tracy's  conterpolnt  criteria  was  satisfactory,  pro- 
vided that  a suitable  ratio  exists  between  the  inner  and  outer  diameters 
of  the  probe.  Since  Bannink  and  Nebbeling  used  a conical  shaped  probe 
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ig.  19.  Cross-Flow  Vector  Velocity  Plot  as  projected  on  F_-  plane  for  Supersonic  Flow  Above 
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in  their  experimental  measurements,  the  shock  location  is  presented  as  a 
band  having  the  width  of  the  observed  pressure  jump.  No  attempt  is  made 
to  define  the  exact  location  of  the  shock  wave  within  this  band.  In 
Figure  21,  it  can  be  seen  that  a close  correlation  exists  between  the 
calculated  and  experimental  locations  of  the  inner  shock  and  sonic  line 
below  £=0.35.  The  small  discrepancy  between  the  measured  and  the  cal- 
culated conical  sonic  line  in  the  upper  region  of  the  flow  field 
(£>■0.35)  is  due  to  the  small  gradients  of  in  the  0 direction,  since 
a small  inaccuracy  in  the  measurements  results  in  a large  r)  variation. 

This  discrepancy  is  indicated  by  the  region  of  uncertainty  shown  in 
Figure  21. 

The  pressure,  temperature,  and  density  contours  in  the  n-£  plane 
are  shown  in  Figures  22-24,  respectively.  These  contour  plots  were 
developed  by  using  a General  Purpose  Contouring  Program  (Ref  134)  on 
the  CDC  6600  computer.  A total  of  780  data  points  were  evaluated  to 
produce  these  figures.  The  internal  shock  wave  and  leading  edge  ex- 
pansion fan  are  clearly  depicted  in  all  three  figures  as  highly  con- 
centrated contour  lines.  In  Figure  22,  the  isobars  indicate  that  the 
inboard  shock  wave,  starting  perpendicularly  from  the  wing  surface, 

» 

extends  into  the  central  region  where  the  expanded  flow  is  dominant. 

In  this  region,  the  internal  shock  weakens  and  eventually  becomes  a 
conical  sonic  line.  Along  the  wing  surface,  the  boundary  layer  is  very 
thin  and  the  pressure  gradient  normal  to  the  surface  is  zero.  Reyond 
the  leading  edge,  a relatively  weak  bow  shock  is  formed  which  diminishes 
very  rapidly  as  it  overflows  the  upper  surface  of  the  wing.  A strong 
temperature  gradient  is  formed  in  the  boundary  layer,  as  seen  in 
Figure  23.  This  gradient  is  slightly  stronger  in  the  region  between 
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the  internal  shock  wave  and  the  leading  edge  expansion  fan.  In  Figure 
24,  it  can  be  seen  that  the  shock  wave-boundary  layer  interaction  is 
very  weak  and  no  flow  separation  occurs  at  the  base  of  the  internal 
shock. 

In  summary,  this  numerical  technique  accurately  predicts  the  basic 
elements  of  the  supersonic  flow  over  the  upper  surface  of  a planar  delta 
wing.  The  numerical  results  compare  quite  favorably  with  Bannink's  ex- 
perimental data,  except  in  the  vicinity  of  the  bow  shock.  These  cal- 
culations show  that  a three-dimensional,  supersonic  viscous  flow  over 
the  expansion  surface  of  a planar  delta  wing  can  be  accurately  approxi- 
mated by  using  a two-dimensional  conical  flow  field  model. 

Hypersonic  Flow.  Cross  (Ref  17)  conducted  a series  of  experimental 
studies  of  the  expansion  side  flow  field  over  a flat  delta  wing  at  hyper- 
sonic speeds.  Several  models  were  used  in  this  investigation.  All  of 
these  models  were  geometrically  similar  flat  plate  delta  wings  with 
sharp  leading  edges  (diameter  approximately  0.003  inch).  The  sweep 
angles  on  these  models  were  73°  and  the  central  wing  chords  varied  from 
3 to  7 inches  with  plate  thicknesses  of  up  to  0.5  inches.  Impact  pres- 
sures were  measured  at  numerous  points  in  the  leeside  flow  field  from 
a=0°  through  01=19°  in  2°  increments.  A 0.040  inch  outside  diameter 
pitot  probe,  mounted  parallel  to  the  upper  wing  surface,  was  used  to 
make  these  pressure  measurements. 

In  order  to  simulate  Cross'  experimental  investigation,  the  fol- 
lowing free  stream  conditions  were  chosen  for  this  calculation 


A constant  stop  size  array,  identical  to  the  upper  half  of  the  grid 
system  in  Figure  9b,  was  used  with  grid  increments  as  shown  in  Table  1. 


Tablel 

Grid  Increments  for  Hypersonic  Flow  Field 
Calcul at  ions 


j 


An 

AJ, 

Grid 

0.019139 

0.00568 

26lix4  5£ 

1 5° 

0.019139 

0.01136 

26nx30£ 

i 

9° 

0.019139 

0.01136 

2bnx30r. 

■ t 

11° 

0.019139 

0.00568 

26r|x76.r, 

f 

15° 

0.019139 

0.00568 

26nx76£ 

The  free  stream  boundaries  (upper  and  right)  were  located  far  enough 
from  the  upper  wing  surface  as  not  to  affect  the  numerical  solutions. 
The  incremental  step  sizes  wer e selected  so  that  no  interpolation  was 
required  In  order  to  compare  experimental  and  numerical  results.  The 
results  of  these  calculations  are  shown  in  Figures  25-35. 

In  Figure  25,  a comparison  is  made  between  t he  experimental  and 
calculated  edge  of  the  viscous  region.  This  viscous  profile  Is 
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determined  by  evaluating  the  Impact  pressures  in  the  ^-direction, 
similar  to  what  is  shown  in  Figure  16.  It  can  be  clearly  seen  (Fig  25) 
that  there  Is  a progressive  increase  in  the  extent  of  the  viscous  re- 
gion as  the  angle  of  attack  is  increased.  At  low  angles  of  attack 
(a<5°)  , the  calculated  profile  can  be  approximated  by  <5~n'-  This  result 
is  similar  to  the  qualitative  flow  behavior  noted  by  Rao  and  Whitehead 
(Ref  131)  in  their  vapor  screen  studies.  However,  from  Cross'  data, 
the  experimental  profile  is  very  nearly  linear.  The  largest  differ- 
ences between  experimental  and  calculated  results  occur  at  the  plane 
of  symmetry,  where  three-dimensional  effects  are  dominant.  For  a>9°,  a 
centerline  "trough"  appears  in  the  calculated  viscous  profile.  This 
trough  is  generated  by  shock-induced  vortices  in  the  viscous  region. 
Cross'  data,  for  a>8°,  shows  a large  region  of  low  impact  pressure  de- 
veloping along  the  centerline  of  the  wing.  This  region  is  located  at 
the  projection  point  of  the  free  stream  velocity  vector  (through  the 
wing  vertex)  on  the  £-r|  plane. 

The  impact  pressure  distribution  for  various  angles  of  attack  and 
^-position  are  shown  in  Figures  26  through  29.  Figure  26  depicts  a 
comparison  between  experimental  and  calculated  impact  pressure  results 
for  et=0°.  The  calculated  values  are  in  good  agreement  with  the  measured 
quantities  except  at  one  point  on  the  upwind  side  of  the  bow  shock  at 
£>'0.0682.  In  this  region,  the  calculated  bow  shock  position  is  slightly 
inboard  of  the  actual  shock  wave.  Figures  27a  and  27b  illustrate  the 
impact  pressure  results  for  a=5  . In  these  figures,  good  agreement  is 
again  shown  between  theory  and  experiment,  except  in  the  vicinity  of  the 
bow  shock  anil  along  the  centerline  ol  the  wing  at  £*0.0682.  This  dis- 
crepancy on  the  symmetry  plane  occurs  because  the  calculated  boundary 
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layer  thickness  Is  slightly  less  than  the  measured  value.  The  same  type 
of  pressure  results  are  also  seen  for  a=9°  (Figs  28a  and  28b),  except 
that  at  5=0. 11505  the  Impact  pressure  disparities  are  much  larger  near 
the  centerline.  These  large  disparities  are  a result  of  the  differences 
between  the  calculated  boundary  layer  trough  and  the  measured  enlarged 
low  pressure  region.  Similar  discrepancies  are  seen  for  a=15°  (Figs 
29a  and  29b)  at  5=0. 15904  and  5=0.20454. 

The  conical  cross-flow  streamlines  for  a=0°,  9°,  and  15°  are  shown 
in  Figures  30  through  32.  At  zero  angle  of  attack  (Fig  30),  a vortical 
singularity  exists  near  the  upper  edge  of  the  viscous  region.  No  flow 
separation  is  seen  along  the  wing  surface  and  no  vortex  is  formed  in 
the  boundary  layer.  As  the  angle  of  attack  is  increased  (as  determined 
in  this  study) , the  vortical  singularity  is  forced  downward  toward  the 
wing  surface  under  the  developing  influence  of  a vortex  in  the  viscous 
region.  At  a=5°,  the  vortical  singularity  is  located  at  the  origin  of 
the  coordinate  system.  All  conical  streamlines  converge  toward  this 
cross-flow  stagnation  point  in  the  r)-5  plane. 

At  a=9°  (Fig  31),  a small  vortex  is  formed  in  the  boundary  layer. 
This  vortex  is  initially  very  close  to  the  wing  surface  and  is  formed 
as  a result  of  shock-induced  boundary  layer  separation.  Because  of  the 
limited  number  of  grid  points  in  the  viscous  region,  only  a coarse  out- 
line can  be  seen  of  this  developing  vortex.  As  the  angle  of  attack  is 
increased,  the  vortex  strength  is  increased  and  the  core  of  the  vortex 
moves  further  above  the  wing.  At  01=11°  and  a=15°  (Fig  32),  the  vortex 
is  well  developed  in  the  viscous  region  and  its  behavior  is  very  similar 
to  the  experimental  observations  of  Cross. 

The  density,  pressure,  and  temperature  contours  in  the  cross-flow 
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Cross-flow  Streamline  Plot  as  Projected  on  F-r;  Plane 
for  Hypersonic  Flow  Above  a Planar  Delta  Wing,  a = 1 


physical  domain  are  shown  in  Figures  3.1  through  38.  These  contour  plots 
depict  the  intricate  flow  field  characteristics  for  a*=0°  to  a*=15°. 
Numerical  results  from  this  study  show  that  several  significant  changes 
occur  in  the  flow  as  the  angle  of  attack  is  increased.  At  n=0°,  the 
cross-flow  is  dominated  by  a strong  lending  edge  shock  wave' due  to 
boundary  layer  displacement  thickness  (Fig  33).  There  Is  no  Internal 
shock  and  thus  the  boundary  layer  remains  attached  on  the  wing  surface. 
The  spanwlse  temperature  distribution,  as  illustrated  in  Figure  34,  in- 
dicates that  the  heat  transfer  is  a minimum  at  the  centerline  and  it 
increases  rather  sharply  at  the  leading  edge.  As  the  angle  of  attack 
is  increased,  the  bow  shock  is  gradually  weakened  by  a developing 
Prandt 1-Meyer  expansion  fan  over  the  leading  edge.  At  0=9°,  an  internal 
shock  is  formed  in  the  inviscid  flow  region.  This  internal  shock  is 
nearly  normal  to  the  wing  surface,  and  at  its  lowest  edge,  is  incident 
upon  the  upper  surface  of  the  viscous  region  (Fig  35).  At  ct=15°,  the 
internal  shock  wave  and  leading  edge  expansion  fan  are  very  strong,  as 
seen  by  the  highly  concentrated  density  contour  lines  in  Figure  36.  The 
bow  shock  is  weak  and  the  pressure  gradient  normal  to  the  wing  surface 
in  the  viscous  region  is  zero  (Fig  37).  Figure  38  shows  a change  in  the 
temperature  profile  from  that  shown  in  Figure  34  for  attached  flow  at 
a=0°.  For  separated  flows,  the  temperature  gradient  (normal  to  the  wing 
surface)  decreases  from  a peak  at  the  centerline,  reaches  a minimum  and 
then  Increases  to  a maximum  at  the  leading  edge.  This  spanwlse  tempera- 
ture behavior  is  similar  to  the  experimental  observations  by  Narayan 
(Kef  25)  at  ot»l5°. 

In  conclusion,  it  can  be  seen  that  the  numerical  method  accurately 
predicts  the  basic  elements  of  leeside  hypersonic  flow  over  a planar 
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Fig.  35  Static  Pressure  Contour  in  Cross 


Flow  Above  a Planar  Delta  Wing. 


delta  wing.  For  the  first  time  in  any  calculation,  the  vortex  develop- 
ment in  the  boundary  layer  of  a delta  wins  is  computed  based  on  a strong 
shock  wave-boundary  layer  Interaction.  The  numerical  results  compare 
quite  favorably  with  Cross'  data  as  well  as  with  the  qualitative  ob- 
servations by  Kao  and  Whitehead  (Kef  131)  and  by  Narayan  (Ref  25).  The 
discrepancies  between  the  calculated  and  experimental  results  are  due 
to  not  modeling  the  compression  side  flow  field  and  not  compensating  for 
the  three-dimensional  effects  around  the  wing  vertex. 


Delta  W 1 ng  Comp res s 1 on  Side  FI oup  Ana 1 v sj s 

The  next  phase  in  this  investigation  is  to  apply  the  conical, 
viscous  flow  approximation  to  a supersonic  flow  over  the  windward  side 
of  a planar  delta  wing.  The  purpose  of  this  effort  is  to  verify  the 
applicability  of  this  method  in  solving  compression  side  flow  fields. 
The  same  simplifying  assumptions,  which  were  used  in  the  expansion  side 
flow  calculations,  are  applied  in  this  calculation. 

For  this  analysis,  only  one  compression  side  flow  field  ease  is 
computed.  The  free  stream  conditions  chosen  for  this  calculation  are 

M - 4.0  A = 50° 

CO 

Re  -=  5.0  x 106  T = 530°K 

X w 


Pr  = 0.72  y *=  1.4 

c*  - 15°  X “ 0.028 


where  T is  equal  to  the  free  stream  stagnation  temperature.  These  flow 

w 

conditions  are  Identical  to  those  used  by  Babaev  (Kef  bb) , Voskresenskli 
(Kef  88),  Ueeman  and  Powers  (Kef  b8) , and  South  and  Klunker  (Kef  b7)  in 
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their  lnvlscld  analyses. 


A 26(r))x45(0 array  was  used  In  this  numerical  calculation.  This 
constant  step  size  array  was  Identical  to  the  lower  half  of  the  grid 
system  shown  In  Figure  9b.  The  n step  size  was  0.059936  with  14  grid 
points  on  the  wing  surface.  The  t,  step  size  was  0.00568.  The  free 
stream  boundary  locations  were  positioned  far  enough  from  the  wing 
surface  and  bow  shock  so  as  not  to  affect  the  numerical  solution.  The 
results  of  this  calculation  are  shown  in  Figures  39-45. 

In  Figure  39,  the  coefficient  of  pressure  on  the  lower  surface  of 
a flat  delta  wing  is  plotted  for  various  spanwise  locations.  The  pri- 
mary area  of  interest,  in  tills  figure,  is  the  subsonic  cross-flow 
region,  where  a variation  in  surface  pressure  occurs.  The  cross-flow 
sonic  line  (as  seen  in  Fig  40)  serves  as  a dividing  line  between  the 
rotational  and  irrotational  portions  of  the  flow.  Close  agreement  is 
seen  between  the  subsonic  numerical  calculations  and  all  of  the  analyti- 
cal solutions,  except  Babaev's  solution.  In  Babaev's  method,  an  attempt 
is  made  to  account  for  the  singularity  which  occurs  at  the  cross-flow 
sonic  point  (on  the  wing  surface).  The  analytical  surface  pressure 
distribution  should  exhibit  a "corner"  or  slope  discontinuity  at  the 
cross-flow  sonic  point.  However,  Babaev's  solution  as  well  as  all  the 
other  analytical  solutions,  show  a very  smooth  pressure  distribution  at 
this  point.  Most  of  the  other  analytical  techniques  ignored  this  weak 
singularity. 

Figures  40  through  43  illustrate  the  cross-flow  Mach  number,  pres- 
sure, temperature,  and  density  contours  in  the  physical  cross-flow  plane, 
respectively.  In  Figure  40,  the  cross-flow  Mach  number  contours  based 
on  the  magnitude  of  the  conical  cross-flow  velocity  components  are  shown. 


Fig.  39.  Spanwise  Pressure  Distribution  on  Compression  Side  of  Planar  Delta  Wing 
at  Supersonic  Speeds 


This  figure  defines  the  cross-flow  subsonic  and  supersonic  regions  and 
the  sonic  line  that  separates  these  regions.  The  sonic  point  is  lo- 
cated at  the  base  of  the  sonic  line  in  the  boundary  layer  or  on  the 


wing  surface  (invlscld  solution  only).  The  shock  wave  from  the  leading 
edge  to  the  sonic  line  is  approximately  planar,  as  seen  in  Figure  41. 

The  invlscld,  analytical  shock  angle,  relative  to  the  leading  edge,  is 
21.3°  compared  with  approximately  1S.51  numerically.  There  is  no 
normal  pressure  gradient  In  the  boundary  layer  and  the  shock  weakens 
slightly  as  it  curves  toward  the  wing  surface  and  approaches  the  plane 
of  symmetry.  In  Figures  42  and  43,  strong  density  and  temperature 

gradients  exist  in  the  thin  boundary  region. 

I 

The  cross-flow  velocity  vectors  and  resulting  conical  streamlines 
are  shown  in  Figures  44  and  45,  respectively.  The  significant  change 
in  magnitude  and  direction  of  the  cross-flow  velocity  vectors  near  the 
shock  wave  indicate  that  the  compression  side  bow  shock  is  very  strong. 
The  streamline  contours  show  that  a vortical  singularity  exists  near  the 
wing  surface  on  the  plane  of  symmetry.  This  is  consistent  with  the 
analytical  results  by  Voskresenski  i (Ref  $$)  and  Melnik  (Ref  Id). 

In  conclusion,  it  can  be  seen  that  this  numerical  technique 
accurately  predicts  the  basic  characterlst ics  of  a compression  side 
flow  field.  The  supersonic  numerical  results  compare  quite  favorably 
with  several  Invlscld  analytical  solutions.  These  calculations  show 
that  this  conical,  viscous  approximat ion  technique  can  be  accurately 
applied  to  high  Reynolds  number,  viscous  compression  side  flow  fields. 

Total  Pelt  a Wing  Flow  Analysis 

The  final  phase  in  this  Invest igat ion  Is  to  examine  and  calculate 
the  total  supersonic  and  hypersonic  flow  field  around  a thin  planar 
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delta  wing.  Previous  calculations  In  this  study  have  demons t rated  the 
applicability  of  the  conical,  viscous  approximation  In  separately  sol- 
ving flows  over  either  the  leeslde  and  or  the  compression  side  of  flat 
delta  wings.  However,  It  Is  still  necessary  to  verify  this  concept  In 
solving  the  total  flow  field  around  a thin  planar  delta  wing.  Thus, 
the  purpose  of  tills  analysis  Is  to  examine  and  evaluate  the  use  ol  the 
conical,  viscous  approximation  In  solving  the  total  delta  wing  flow 
field.  In  this  calculation,  the  upper  and  lower  surface  flow  fields 
are  allowed  to  interact  with  each  other  in  reaching  a steady  state  solu- 
tion. However,  the  three-dimensional  effects  generated  by  the  wing 
vertex  are  still  neglected  in  solving  this  flow  problem. 


Supersonic  Flow.  For  the  supersonic  flow  analysis,  t lie  experimen- 
tal test  conditions  used  by  Hannlnk  and  Nebbeling  (Ret  1 r0  are  applied 
in  this  calculation.  These  flow  conditions  are 
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where  the  Reynolds  number  is  based  on  a root  chord  length  ot  0.171  tt  . 
This  Is  the  same  case  examined  tor  the  leeside-onlv  solution. 

The  cornput at i ona 1 domain,  as  shown  In  Figure  <*b , consisted  ot  a 
26(n)xV>(f.)  grid  array  with  ?'!  rows  of  grid  points  above  and  below 


1?8 


A 


the  wing  surface.  The  r)  step  size  was  0.063158  with  14  grid  points  used 
to  model  the  wing.  A constant,  but  different,  £ step  size  was  used  for 
each  measurement  height  (pitot  pressure  measurements  by  Ranntnk,  Ref  15) 
so  that  experimental  and  numerical  results  could  be  compared  without 
Interpolation.  The  free  stream  boundaries  were  located  far  enough  from 
the  wing  surface  so  as  not  to  affect  the  numerical  solution.  The  re- 
sults of  these  calculations  are  shown  In  Figures  46-52. 

In  Figures  4b  and  47,  several  spanwlse  pitot  pressure  measurements 
are  shown  for  various  heights  above  the  delta  wing.  The  numerical  re- 
sults are  compared  with  the  Invlscld  solutions  by  Kut ler  (Ref  89)  and 
with  the  experimental  data  by  Bannlnk  and  Nebbeling  (Ref  15).  It  can 
be  seen  that  the  total  numerical  flow  field  solution  Is  similar  to  that 
of  the  expansion-side-only  supersonic  solution.  The  only  significant 
difference  Is  that. the  bow  shock  In  the  leeslde  flow  field  Is  slightly 
stronger  for  £ < 0.3626.  This  Is  due  to  the  influence  of  the  lower 
surface  shock  wave  as  It  encircles  the  leading  edge  of  the  wing.  The 
coarse  grid  spacing  and  the  use  of  the  numerical  shock  capturing  tech- 
nique results  in  smearing  the  shock  over  several  grid  points.  This 
causes  the  leeslde  shock  to  be  displaced  slightly  upstream  of  the  meas- 
ured bow  shock  location.  A similar  error  was  noted  by  Bazzhln  (Ref  92) 
in  his  delta  wing  calculations. 

Figure  48  depicts  the  conical  cross-flow  velocity  components  pro- 
jected on  the  n-£  plane.  On  the  compression  side,  the  bow  shock  is 
clearly  delineated  by  an  abrupt  change  in  strength  and  direction  of  the 
velocity  vectors.  The  locus  of  the  velocity  vectors  trace  out  "pseudo" 
streamlines  which  converge  at  two  vortical  singularity  points.  These 
vortical  singularity  points  are  located  above  and  below  the  wing 
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Numerical  Result 


Fig.  46.  Pitot  Pressure.  Distribution  in  Spanwise  Direction  for 
Supersonic  Flow  Above  a Planar  Delta  Wing,  C “ 0.0718. 
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Fig.  48.  Cross-flow  Velocity  Vector  Plot  as  Projected  on  \-r  Pl^ne 
for  Supersonic  Flow  Around  a Planar  Delta  Wing,  cx  = 12 


surface  near  Che  origin  of  the  coordinate  system.  Figure  49  illustrates 
some  of  the  conical  cross-flow  streamlines  (1550  points)  and  how  these 
steady  state  path  lines  converge  toward  the  two  vortical  singularity 
points.  These  streamline  results  are  consistent  with  the  inviscid, 
analytical  solutions  by  Melnik  (Ref  19)  and  with  the  experimental  ob- 
servations by  Bannlnk  and  Nebbeling  (Ref  15). 

The  pressure,  temperature,  and  density  contour  plots  for  the  physi- 
cal cross-flow  plane  are  shown  in  Figures  50-52.  The  wing  thickness  is 
exaggerated  in  these  figures  in  order  to  clearly  depict  its  location.  A 
total  of  1550  data  points  were  used  to  produce  these  contours  including 
a set  of  double-values  grid  points  representing  the  wing  surface.  In 
Figure  50,  the  isobars  indicate  that  a strong  shock  wave  is  formed  below 
the  delta  wing.  This  shock  wave  weakens  slightly  as  it  bends  toward  the 

wing  surface  and  approaches  the  plane  of  symmetry.  The  influence  of  this 

compression  wave  is  felt  above  the  delta  wing,  as  seen  by  the  contour 

lines  encircling  the  leading  edge  of  the  wing.  A large  leading  edge  ex- 

pansion fan  is  formed  above  the  wing  surface  and  inboard  of  the  bow  shock. 
This  expansion  fan  weakens  the  leeside  internal  shock  wave  and  eventually 
transforms  it  into  a conical  sonic  line.  Along  the  wing  surfaces,  the 
boundary  layers  are  very  thin  and  the  pressure  gradients  normal  to  these 
surfaces  (in  the  £j  direction)  are  zero. 

As  seen  in  Figure  51,  a strong  temperature  gradient  exists  in  both 
boundary  layers.  This  temperature  gradient  (normal  to  the  wing  surface) 
is  particularly  strong  near  the  leading  edge.  In  Figure  52,  it  can  be 
seen  that  the  leeside  internal  shock  wave-boundary  layer  interaction  is 
very  weak  and  that  no  flow  separation  occurs  at  the  base  of  the  shock 
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rig.  52.  Static  Density  Contour  in  Cross-flow  Plane  for  Supersonic  Flow 
Around  a Planar  Delta  Wing. 


In  summary,  Che  numerical  solution  of  the  supersonic  flow  around  a 
thin,  planar  delta  wing  provides  a more  accurate  and  complete  solution 
of  the  delta  wing  flow  field  than  does  the  leeside-only  solution.  This 
total  flow  field  solution  contains  all  the  basic  elements  of  the  flow 
and  these  results  compare  quite  favorably  with  Bannink's  experimental 
data.  The  discrepancies  noted  in  comparing  numerical  and  experimental 
results  are  due  to  three-dimensional  effects  in  the  viscous  region,  use 
of  a coarse  grid  in  the  computational  plane,  and  not  properly  modeling 
the  delta  wing  thickness.  However,  the  results  indicate  that  the  three- 
dimensional  flow  field  can  be  approximated  by  using  a two-dimensional 
conical,  viscous  flow  field  model  which  would  be  useful  in  design 
applications . 

Hypersonic  Flow.  For  the  hypersonic  flow  field  analysis,  the  free 
stream  and  surface  boundary  conditions  chosen  for  this  calculation  are 


M = 10.17 

T 

= 1780°R 

oo 

o 

Re  = 3.345  x 105 

P 

= 596  psia 

X 

o 

Pr  = 0.72 

y = 

1.4 

X = 1-8 

T 

w 

= 1259. 7°R 

a = 0°,  5°,  9°.  11°,  and  15° 

A = 

75° 

where  the  Reynolds  number  is  based  on  a reference  length  of  x = 5.5  sec  a 
inches.  These  flow  conditions  are  the  same  as  those  used  by  Cross 
(Ref  17)  in  his  experimental  studies  of  the  hypersonic  flow  over  the 
lecside  of  a flat  delta  wing. 

In  this  investigation,  a constant-step  size  array  (Fig  9b)  was  used 
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to  calculate  the  flow  field.  This  computational  domain  consisted  of  a 
26(r|)xl20(O  grid  array  with  75  rows  of  £ grid  points  above  the  wing  and 
44  rows  below  the  wing.  The  n step  size  was  0.019139  with  14  grid  points 
on  the  wing  surface.  The  £ step  size,  for  all  the  hypersonic  cases,  was 
0.00568.  The  exterior  or  free  stream  boundaries  were  located  far  enough 
from  the  wing  surface  so  as  not  to  affect  the  numerical  solution.  The 
results  of  these  calculations  are  shown  in  Figures  53-73. 

In  Figure  53,  a comparison  is  made  between  the  measured  and  cal- 
culated edge  of  the  viscous  region.  This  boundary  layer  profile  is  de- 
termined by  evaluating  the  impact  pressure  distribution  in  the  in- 
direction similar  to  what  is  done  in  Ref  17.  For  a < 5°,  the  edge  of 
the  calculated  viscous  region  is  very  similar  to  that  predicted  by  two- 
dimensional  laminar  boundary  layer  theory.  The  computed  boundary  layer 
thickness,  5,  for  a = 0°,  is  identical  to  that  calculated  for  the  ex- 
pansion-side-only solution;  however,  for  a = 5°,  6 is  greater  and  more 
accurate  than  the  leeside-only  result.  At  a = 9°,  a centerline  trough 
is  computed  in  the  cross-flow  viscous  profile.  This  trough  is  generated 
by  a developing  shock-induced  vortex  in  the  boundary  layer.  For  a 11°, 
a large  viscous  bubble  or  "hump"  appears  in  the  plane  of  symmetry  of  the 
boundary  layer  edge.  This  viscous  bubble  occurs  as  a result  of  shock- 
induced  separation  behind  a strong  leeside  bow  shock.  This  strong  bow 
shock  is  clearly  depicted  in  the  contour  plots  shown  later  in  Figures 
71  and  72.  Good  qualitative  agreement  is  seen  between  the  experimental 
and  calculated  results  even  though  the  computed  boundary  layer  thick- 
ness is  slightly  greater  than  the  measured  value. 

In  Figure  54,  the  bow  shock  location  on  the  expansion  side  of  the 
delta  wing  is  shown  for  various  angles  of  attack.  For  n ^ 0.1,  the 


139 


calculated  shock  profile  Is  slightly  less  than  the  measured  shock  loca- 
tion. This  discrepancy  occurs  because  the  displacement  thickness  of  the 
thin,  planar  delta  wing  Is  smaller  than  that  of  the  actual  wind  tunnel 
model.  Along  the  plane  of  symmetry,  good  agreement  Is  seen  between  the 
calculated  and  measured  shock  shape  for  ot  £ 5l  . However,  for  a > 11°, 
the  numerical  shock  profile  Is  slightly  greater  than  the  experimental 
value.  This  error  is  due  to  the  displacement  effects  of  the  cal- 
culated boundary  layer  In  the  symmetry  plane. 

The  impact  pressure  distribution  for  various  angles  of  attack  and 
^-positions  Is  shown  in  Figures  55  through  58.  Figure  55  illustrates 
the  comparison  between  experimental  and  calculated  impact  pressure  re- 
sults for  a = 0°.  These  results  are  similar,  but  more  accurate  than  the 
leeside-only  results,  particularly  at  £=0.1l3b.  Along  the  plane  of 
symmetry,  the  calculated  boundary  layer  thickness  is  slightly  less  than 
the  measured  value.  This  discrepancy  is  very  apparent  at  £=0.0682, 
where  a large  difference  exists  between  the  measured  and  calculated 
impact  pressures.  Figures  56a  and  56b  depict  the  impact  pressure  re- 
sults for  ct  » 5°.  In  these  figures,  a closer  agreement  is  again  seen 
between  theory  and  experiment.  The  calculated  pressure  ratio  across  the 
bow  shock  Is  stronger  than  that  of  the  leeside-only  solution,  but  weaker 
than  the  experimental  results.  The  flow  resolution  around  the  shock 
wave  is  poor  due  to  shock  smearing  In  the  coarse  computational  grid. 
Similar  results  are  seen  for  it  = 9*  (Figs  57a  and  57b)  and  for  a - IS1 
(Figs  58a  and  58b).  However,  for  a = 15C  , the  three-dimensional 
effects  near  the  centerline  are  more  dominant,  as  seen  at  (,*0.15904  and 
(,-0. 20454. 

The  conical  cross-flow  velocity  vectors  and  streamline  plots  tor 


142 


0.02301 


0.11505 

0.20709 


A 


A 


A 


/ 


/ 


/ 

/ 

-O-0-"  n 


a ■ 0°,  9°,  11°,  and  15°  are  shown  in  Figures  59  through  65.  At  zero 
angle  of  attack  (Figs  59-60) , vortical  singularities  occur  near  the 
edges  of  the  upper  and  lower  viscous  regions.  No  flow  separation  is 
seen  along  the  wing  surface  and  no  vortex  is  formed  in  the  boundary 
layer.  At  a - 5°  (as  determined  in  this  study),  these  vortical  singular- 
ities are  forced  toward  the  wing  surface  on  both  sides  of  the  wing.  A 
further  increase  in  angle  of  attack  to  a * 9°,  results  in  a detached 
vortical  singularity  below  the  wing  surface  (Fig  61) . The  flow,  on 
both  sides  of  the  wing,  remains  attached  although  there  is  an  increase 
in  circulation  (the  beginnings  of  a vortex)  on  the  leeside  near  the 
origin  at  a = 9°. 

At  a = 11°,  (Figs  62-63),  a large  reversed  flow  region  or  vortex  is 
formed  in  the  boundary  layer.  This  vortex  is  clearly  depicted  in  the 
conical  streamline  plot  in  Figure  63.  Above  this  vortex  near  the  edge 
of  the  leeside  boundary  layer  is  a vortical  singularity  point.  Another 
vortical  singularity  point  is  formed  below  the  wing  surface  near  the  edge 
of  the  viscous  region.  The  leeside  cross-flow  separation  point  is  at 
ri-0.1206  with  the  experimental  point  at  n=0.1166.  As  the  angle  of  attack 
is  increased,  the  vortex  strength  is  also  increased  and  the  core  of  the 
vortex  moves  further  above  the  wing.  At  a = 15°  (Fig  64)  the  cross-flow 
separation  point  is  at  r)=0.1301  while  the  experimental  separation  point 
is  at  ri=0.1432.  The  leeside  vortical  singularity  moves  further  from  the 
wing  surface  as  the  boundary  layer  thickness  increases  (Fig  65).  This 
behavior  of  the  vortical  singularity  points  is  similar  to  that  seen  in 
Tracy's  experimental  studies  on  a cone. 

The  pressure,  density,  and  temperature  contours  in  the  cross-flow 
physical  plane  are  shown  in  Figures  66  through  73.  These  figures 
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Fig.  59.  Cross-flow  Velocity  Vector  Plot  ns  Pro|ecLeil  on  f.-n  Plane 
for  Hypersonic  Flow  Around  n I’lnnnr  Delta  King,  Oi  " 0 . 
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Fig.  61.  Cross-flow  Velocity  Vector  Plot  as  Projected  on  f’-n  P^ane 
for  Hypersonic  Flow  Around  a Planar  Pelta  Wing,  a = 9 . 
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Fip.  64.  Cross-Flow  Velocity  Vector  Plot  as  Projected  on  ’)  Pl^ne 
for  Hypersonic  Flow  Aroutul  a Planar  Pol  to  Win);,  <x  “ 15  . 
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Cross-flow  Streamline  Plot  as  Projected  on  f,-n  Plain 
for  Hypersonic  Flow  Around  a Planar  Delta  Wlnp,  ct  ■■ 
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illustrate  the  characteristics  of  the  total  flow  field  for  a = 0°  to 
a ■ 15°.  At  zero  angle  of  attack,  the  cross-flow  is  dominated  by  a 
strong  leading  edge  shock  wave  (Fig  66) . This  shock  wave  is  formed 
because  of  the  displacement  effects  of  the  boundary  layer.  The  spanwise 
temperature  distribution  (Fig  67)  in  the  boundary  layer  indicates  that 
the  heat  transfer  rate  gradually  increases  from  a minimum  at  the  center- 
line  to  a maximum  at  the  leading  edge  of  the  wing.  There  is  no  internal 
shock  wave  and  thus  the  boundary  layer  remains  attached  to  the  wing  sur- 
face. As  the  angle  of  attack  is  increased,  the  pressure  ratio  across 
the  compression-side  shock  wave  also  increases.  This  shock  wave  has  a 
strong  influence  on  the  leeside  bow  shock,  as  seen  in  Figures  68-69. 

At  a = 9°,  an  internal  shock  is  formed  on  the  leeside  of  the  delta  wing. 
This  internal  shock  is  normal  to  the  wing  surface  and,  at  its  lowest 
point,  is  incident  upon  the  upper  edge  of  the  viscous  region  (Fig  70). 

At  a = 11°  (Fig  71) , the  interaction  of  the  internal  shock  wave  and  the 
boundary  layer  is  so  strong  that  a vortex  is  formed  in  the  viscous  re- 
gion. The  leeside  bow  shock  and  Prandtl-Meyer  expansion  fan  are  much 
stronger  than  those  calculated  in  the  leeside-only  solution.  At  a = 15°, 

i 

the  bow  shock  and  expansion  fan  are  the  dominant  features  of  the  upper 
surface  flow  field  (Fig  72) . Although  an  internal  shock  wave  does  exist, 
this  shock  is  much  weaker  than  the  internal  shock  calculated  for  the 
leeside-only  solution.  The  pressure  gradient  (normal  to  the  wing  sur- 
face) in  the  boundary  layer  is  zero,  except  in  the  separated  flow  re- 
gions.  The  spanwise  temperature  distribution  (Fig  73),  for  the  separa- 
ted flow,  is  similar  to  that  computed  for  the  expansion-side-only 
solution.  The  only  major  difference  is  that  the  heat  transfer  rate  for 
the  total  flow  field  is  slightly  lens  near  the  leading  edge. 

158 


I 


I 


161 


Fig.  69.  Static  Pressure  Contour  In  Cross-flow  Plane  for  Hypersonic 
Flow  Around  a Planar  Delta  Wing,  at  = 9°. 
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Fig.  71.  Static  Pressure  Contour  In  Cross-flow  I’l.inc  for  Hypersonic 
Flow  Around  a Planar  Delta  Wing,  a = 11  . 


In  conclusion,  it  can  be  seen  that  the  total  flow  field  solution  is 
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more  accurate  and  complete  than  the  leeside-only  solution.  These 
numerical  results  indicate  that  the  compression  side  flow  field  has  a 
significant  effect  on  the  flow  characteristics  on  the  leeside  of  the 
wing.  The  pressure  ratio  across  the  upper  surface  bow  shock  is  stronger 
than  that  found  in  the  leeside-only  solution.  This  results  in  a larger 
region  of  reversed  flow,  a more  accurate  prediction  of  the  cross-flow 
separation  point,  and  the  calculation  of  a viscous  bubble.  The  hyper- 
sonic numerical  results  compare  quite  favorably  with  Cross'  data  as  well 
as  with  the  qualitative  observations  by  Rao  and  Whitehead  (Ref  131)  and 
by  Narayan  (Ref  25) . The  discrepancies  between  the  calculated  and  meas- 
ured results  are  due  to  three-dimensional  effects  in  the  viscous  region 
and  not  properly  modeling  the  delta  wing  thickness.  However,  these  re- 
sults indicate  that  the  conical  viscous  approximation  technique  can  be 
used  to  adequately  predict  the  three-dimensional  flow  around  a thin  delta 
wing. 

Computational  Statistics 

The  two  most  important  statistics  for  any  computer  program  are  the 

execution  time  per  point  and  per  time  step  and  the  storage  requirements. 

These  statistics  are  used  to  measure  the  efficiency  of  the  numerical 

codes  and  the  size  of  computer  needed  to  run  these  programs.  In  this 

investigation,  the  computer  programs  were  run  on  the  Control  Data 

Corporation  6600  computer.  The  execution  time  for  the  various  numerical 

-3  -3 

solutions  ranged  from  3.07  x 10  to  3.73  x 10  sec/grid  point/time 
step.  The  computer  storage  requirements  and  the  numerical  damping 
constants  varied  for  each  case.  These  values  are  shown  in  Table  II. 
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Table  II 


Computer  Storage  Requirements  and 
Numerical  Damping  Constants 


Case 

Grid 

Storage 
(K  words) 

Damping 

Constants 

Cone 

63rix30£ 

167. 3K 

No  Damping 

Supersonic 
Upper,  Delta 
Wing 

26r|x30£ 

101K 

No  Damping 

Hypersonic 
Upper,  Delta 
Wing 

26r|x45£  (ct=0°) 

26nx30£  (ct=5°.  9 ) 
26rjx76£  (a=ll  , 15°) 

125. 2K 

102. 2K 

174. 4K 

3=1.0 

cj=0.06 

Cj  =0 . 0 

Supersonic 
Lower,  Delta 
Wing 

26nx45£ 

125. IK 

6=1.0 

cA=0.12 

Cj=0.12 

Supersonic 
Total,  Delta 

26nx59£ 

151. 5K 

3=1.0 

ci=0.06 

cj=0.0 

Hypersonic 
Total,  Delta 
Wing 

26nxl20£ 

167. 3K 

6i=6,=l-0 
83=  -20.0 
cj=0.40 

Cj=0. 10 


The  total  execution  time  for  each  case  varied  from  3 to  12  hours, 
depending  on  how  the  stability  criteria  (At)  was  used,  the  use  and 
magnitude  of  the  numerical  damping  terras,  the  size  of  the  computational 
mesh,  the  convergence  criteria,  and  the  characteristics  of  the  physical 
flow  field. 
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VI.  Conclusions  and  Recommendations 


A numerical  method  was  used  to  compute  the  supersonic  and  hyper- 
sonic, viscous  flow  fields  around  a thin,  planar  delta  wing.  These 
solutions  were  obtained  by  solving  the  unsteady  governing  equations 
subject  to  a conical  approximation.  The  integration  technique  used  was 
the  second-order  accurate  MacCormack  finite-difference  scheme.  This 
integration  was  performed  on  a constant  step  size  array  generated  by  a 
conical  coordinate  transformation.  The  solutions  obtained  were  for  a 
Mach  number  range  of  2.94  to  10.17,  a local  Reynolds  number  range  of 
3.345  x 10^  to  5.0  x 10^,  and  angles  of  attack  from  -15°  to  + 15°. 
Numerical  oscillations  in  these  solutions  (as  a result  of  shock 
capturing)  were  reduced  by  applying  normal  stress  damping  and  a fourth- 
order  density  damping  term  to  the  finite-difference  equations.  A 
stability  criteria  (maximum  At)  was  computed  and  used  based  on  an 
analysis  of  the  linearized  governing  equations.  The  numerical  results 
were  compared  with  experimental  data  (Ref  15-17),  various  analytical 
solutions  (Ref  66-68,  88,  and  89),  and  several  qualitative  observations, 
such  as  vapor  screen  and  oil  flow  techniques  (Refs  25  and  131).  From 
these  results,  the  following  significant  conclusions  were  drawn  based 
on  the  present  investigation: 

(1)  This  numerical  technique  accurately  predicts  the  supersonic 
flow  around  a thin,  planar  delta  wing,  with  supersonic  leading  edges. 
Good  agreement  was  obtained  between  calculated  results  and  the  ex- 
perimental data  by  Bannink  (Ref  15).  All  the  basic  elements  of  the 
flow  field  (i.e.  shock  waves,  boundary  layers,  and  sonic  lines)  wore 
calculated  and  were  found  to  be  essentially  correct  in  magnitude  and 


location.  For  the  first  time  In  any  reference,  the  shock  wave-boundary 
layer  Interaction  was  computed  and  Its  effects  seen  In  the  Impact  pres- 
sure profile  for  £ ■ 0.0718  (Figures  17  and  46).  The  numerical  solutions 
for  the  leeside-only  and  the  total  flow  fields  were  found  to  be  almost 
Identical,  except  near  the  bow  shock.  In  this  region,  the  difference  In 
calculated  results  was  due  to  not  modeling  the  compression  side  flow 
field  for  the  upper-surface-only  solution.  The  largest  discrepancies 
between  calculated  and  measured  impact  pressures  also  occurred  near  the 
leeside  bow  shock.  The  maximum  impact  pressure  error  for  the  leeside- 
only  solution  was  12.0%  while  that  for  the  total  flow  field  was  4.0%. 
These  discrepancies  were  attributed  to  use  of  a coarse  grid  in  the  com- 
putational domain,  improper  modeling  of  the  delta  wing  thickness,  and 
neglecting  the  lower  surface  flow  field  influence  (leeside-only  solu- 
tion). From  these  numerical  results,  it  can  be  seen  that  the  three- 
dimensional  supersonic  flow  field  can  be  accurately  approximated  by 
using  a conical,  viscous  flow  field  model. 

(2)  This  numerical  method  adequately  modeled  the  inviscid  and 
viscous  hypersonic  flow  around  a thin  delta  wing.  Solutions  were  ob- 
tained for  both  the  upper  and  total  flow  fields.  These  solutions  com- 
pared quite  favorably  with  Cross'  data  as  well  as  with  the  qualitative 
observations  by  Rao  and  Whitehead  (Ref  131)  and  by  Narayan  (Ref  25) . 

The  shock  wave-boundary  layer  interaction  was  accurately  approximated 
and,  for  the  first  time,  an  embedded  vortex  was  computed  in  the  viscous 
region.  A significant  difference  was  seen  between  the  calculated 
leeside-only  and  total  flow  field  solutions  (see  Figs  32  and  65).  This 
difference  was  due  to  the  strong  interaction  which  occurs  between  the 
upper  and  lower  flow  fields.  In  the  total  flow  field  solution,  the 
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detached  bow  shock  and  leading  edge  expansion  fan  were  found  to  be  the 
dominant  features  In  the  upper  surface  flow  field.  These  flow  character- 
istics caused  a weak  internal  shock  wave  and  a large  reversed  flow  re- 
gion to  form  above  the  wing  at  a > 11°.  For  the  first  time  In  any  cal- 
culation, the  cross-flow  separation  point  was  accurately  predicted 
(within  3. OX)  and  the  viscous  "bubble"  in  the  symmetry  plane  was  com- 
puted. The  total  flow  field  solution  provided  an  improved  insight  into 
the  behavior  of  the  cross-flow  vortical  singularities  and  a more  accurate 
description  of  the  flow  field  than  the  leeside-only  solution. 

(3)  A very  effective  technique  was  used  to  calculate  the  surface 
conditions  at  the  leading  edge  of  the  total  delta  wing.  This  method 
permitted  the  flow  properties  to  be  triple  valued  at  the  singularity 
point.  The  standard  boundary  conditions  of  non-slip  at  the  surface  and 
an  isothermal  surface  were  used.  The  normal  momentum  equations  in  the 
£ and  n directions  were  used  to  calculate  the  upper,  side,  and  lower 
values  of  pressure.  The  three  values  of  density  were  determined  from 
the  equation  of  state.  This  numerical  modeling  of  the  leading  edge  or 
singularity  point  produced  a stable  and  accurate  solution  of  the  flow 
field  in  the  vicinity  of  this  nodal  point. 

(A)  A satisfactory  stability  criteria  analysis  was  performed  on 
the  finite-difference  form  of  the  linearized  governing  equations.  This 
analysis  accounted  for  both  the  inviscld  and  viscous  dominant  regions 
in  the  flow  field.  A maximum  allowable  time  step  was  determined  and 
used  in  the  numerical  Integration  of  the  governing  equations. 

(5)  A normal  stress  damping  term  and  a fourth-order  density  damp- 
ing term  were  successfully  incorporated  into  the  numerical  integration 
procedure.  The  normal  stress  damping  was  used  to  control  the  initial 


transients  due  to  ill-suited  Initial  conditions.  The  density  dumping 


term  was  used  to  reduce  numerical  oscillations  around  shock  waves  and 
expansion  fans.  The  magnitude  of  both  damping  terms  was  set  so  as  to 
be  effective  in  the  Invlscld  flow  region  but  not  to  change  the  viscous 
region  behavior  or  modify  tiie  effective  Reynolds  number,  appreciably. 

Thus,  it  can  be  seen  that  tills  numerical  technique  accurately  pre- 
dicts the  supersonic  and  hypersonic  flow  fields  around  a thin  planar 
delta  wing.  These  numerical  results  compare  quite  favorably  with  ex- 
perimental data,  various  analytical  solutions,  and  several  qualitative 
observations.  This  Investigation  has  demonstrated  the  feasibility  of 
using  the  conical  flow  approximation  in  calculating  the  vlscous- 
inviscld  flow  fields  around  thin  delta  wings  for  x < 0(1). 

Further  research  in  this  area  is  still  needed  in  order  to  examine 
the  total  spectrum  of  supersonic  and  hypersonic  flows  around  a thin 
delta  wing.  Several  recommendations,  bv  the  author,  are  proposed  for 
future  work.  These  include  the  following: 

(1)  A faster  and  more  efficient  numerical  algorithm  should  be  used 
to  solve  this  flow  field  problem.  Sluing  (Ret  l U>)  lias  developed  a time- 
dependent  implicit-explicit  hybrid  scheme  which  is  S.b  times  faster  than 
the  current  explicit  method.  This  hybrid  scheme  could  be  incorporated 
into  the  current  computer  codes  in  order  to  produce  a faster,  but 
equivalent  (accuracy  within  4%),  result. 

(2)  A conical,  body-fitted,  curvilinear  coordinate  system  (Ref 
137)  should  be  used  to  solve  i lows  around  thin  delta  wings  of  various- 
conical  cross-sect  ions . In  this  coordinate  system,  the  grid  system  is 
generated  by  solving  a set  of  elliptical  partial  dltferential  equations, 
with  one  coordinate  being  coincident  with  each  boundary  contour  in  the 
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physical  domain.  Use  of  this  type  of  coordinate  transformat  Ion  would 


result  In  solutions  about  delta  wings  of  various  cross-sections  with 
the  same  degree  of  ease  as  for  planar  delta  wings. 

(3)  As  a first  attempt  at  calculating  the  turbulent  flow  around 
a planar  delta  wing,  an  eddy  viscosity  model  should  be  Incorporated 
Into  the  governing  equations.  Several  algebraic  turbulent  models  could 
be  used  in  this  study,  provided  the  rate  of  strain  Is  fairly  small  and 
the  deviation  of  the  principal  axes  between  the  Reynolds  stress  tensor 
and  the  strain-rate  tensor  is  small  (Ref  138).  This  effort  should  be 
attempted  only  whet,  an  adequate  turbulent  experimental  data  base  exists. 

(4)  Finally,  additional  experimental  data  is  needed  on  supersonic 
and  hypersonic  flow  fields  around  delta  wings.  A more  detailed  ex- 
amination is  required  of  the  viscous  region  at  various  angles  of 
attack,  and  Reynolds  numbers.  An  experimental  study  should  be  pursued 
to  investigate  the  limits  of  the  conical  flow  approxlmat ion  In  evalua- 
ting flows  around  delta  wings.  More  experimental  measurements  are 
needed  in  order  to  examine  the  behavior  of  vortical  singularities  in 
the  cross-flow  of  a low  Reynolds  number,  hypersonic  flow  field. 
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APPENDIX  A 

Trans format  ton  Derivatives 


This  appendix  contains  a comprehensive  set  of  transformation 
derivatives  used  In  the  numerical  integration  of  the  governing  equations. 
Since  the  intent  here  Is  to  provide  a quick  reference  only,  most  of  the 
algebraic  development  is  omitted. 

Transformation  Derivatives  for  the  Cone 


For  the  flow  field  calculations  over  a cone,  the  coordinate  trans- 
formation is 


Transformation  Derivatives  for  the  Delta  Winn 

For  the  flow  around  a dolta  wing,  the  coordinate  transformation  Is 


The  transf ormat Ion  derivatives  are 
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APPENDIX  B 

Convergence  of  Iterative  Mot  hods 
and  Dotormln.it  ton  of  It  ora  t Ion 
Errors 

All  iterative  techniques  generate  a sequence  of  numbers  or  vectors. 

These  iterative  methods  are  convergent  If  their  sequences  converge  to 

the  numbers  or  vectors  which  satisfy  the  given  problem.  Therefore, 

given  an  initial  vector  X*'0^,  where  X = Jxj,  x.,,  . . . x^J.an  iterative 

(k) 

technique  generates  a sequence  of  vectors  X which,  hopefully,  con- 
verges to  a limit  vector,  A. 

A fundamental  theorem  of  numerical  analysis  asserts  that  a sequence 

(k) 

converges  If  and  only  if  it  Is  a Cauchy  sequence.  A sequence  X is 
a Cauchy  sequence. if  for  every  c '•0  there  exists  a positive  number  N 
such  that  for  all  integers  n,  m >N  we  have 

| |x(«)_  | | < e (Bl) 

where  ||  | | is  some  vector  norm.  The  choice  of  the  particular 

norm  is  not  important  since  it  can  be  shown  that  all  vector  norms  are 
equivalent  (Ref  143). 

The  importance  of  the  Cauchy  property  is  that  the  convergence  of  a 
sequence  can  be  ascertained  without  knowledge  of  the  limit  vector,  A. 

The  significance  of  this  concept  with  regards  to  numerical  methods  in 
which  A is  not  known  is  readily  apparent. 

In  this  Investigation,  the  numerical  convergence  method  is  derived 
from  the  definition  of  a Cauchy  sequence  and  is  referred  to  as  the 
Cauchy  method.  Numerical  iterations  are  carried  out  until  a condition 
of  the  form 
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4>  -4' 


"ij  . - 

V t 
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la  mot  where  4>  - u,  v.  w.  p,  p.  o and  e Is  chosen  positive  number.  4^ 

is  defined  as  4>|k]  for  all  the  primitive  variables,  except  for  the 

veloc It v terms  where  V Is  used.  The  local  error,  E(x.  , x,,  . . .x  k)  , 
• max  1 •-  11  • 


at  iterat Ion  k Is  defined  as  the  vector  norm 
the  true  error  at  each  Iteration  k Is 


(k'  v(k)  . 

• *i.j  ' A 


L(k)_  x(k-n 


wh  l le 


t.J 


(113) 


The  number  t appearing  In  l'q  it.'  is  called  the  convergence  criteria. 

The  condition  specified  by  Kq  It.’  Is  that  t or  the  Iterative  convergence 
of  the  finite-difference  technique.  The  value  of  the  convergence 
criteria  for  this  Investigation  is  10  s. 


APPENDIX  c: 

Stab  1 1 1 tv  Analysis 

The  theoretical  Investigation  of  stability  is  particularly  complex 
for  the  difference  schemes  associated  with  the  governing  equations. 
Approaches  to  studying  this  phenomenon  are  reported  In  Koache  (Kef  117) 
and  Rlchtmyer  and  Morton  (Ref  109).  An  approximate  method  which  tends 
to  yield  the  best  results  for  a set  of  general  nonlinear  equations  is 
the  amplification  matrix  theory  by  von  Neumann.  This  method  consists 
of  examining  the  linearised  difference  equations  for  the  amplification 
of  short  wavelength  oscillations  superimposed  on  an  exact  solution. 

The  growth  and  decay  of  these  oscillations  can  be  predicted  by  applying 
the  harmonic  analysis  of  von  Neumann.  In  his  analysis,  the  boundary 
conditions  have  no  effect  on  the  stability  result  and  the  exact  solu- 
tion of  the  governing  equations  is  smooth.  This  latter  assumption 
allows  the  coefficients  of  the  partial  differential  equations  to  be 
treated  as  constants  (locally).  The  stability  conditions  predicted  bv 
this  theory  result  in  a local  stability  condition  which  places  a bound 
on  the  time  Increment  used  In  the  numerical  integration. 

Richtmver  and  Morton  (Ref  109)  showed  that  the  von  Neumann  stabil- 
ity conditions  for  the  nonconservative  form  of  the  governing  equations 
is  the  same  as  that  for  the  conservative  form.  Since  the  analysis  is 
easier  for  the  nonconservat Ive  form,  the  following  nond tmens iona l i red 
governing  equations,  as  suggested  bv  MacCormack  and  Raldwin  (Ret  ITT), 
were  considered: 
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elements  of  these  matrices 

can  be 

found 

in  App 

end  lx 

E. 

In  the  above  equations,  the  Prandtl  number,  the  Reynolds  number, 
and  the  nodal  point  locations  are  assumed  to  be  constant.  In  addition, 
the  dissipation  terms  of  the  energy  equation  which  contain  quadratics 
in  first-order  derivatives  have  been  deleted  based  on  an  analysis  by 
Kentzer  (Ref  124) . 

The  stability  criteria  for  the  governing  equations  are  determined 
by  examining  three  distinct  parts  of  these  equations:  the  inviscld, 
diffusion,  and  mixed-derivative  parts.  The  magnitude  of  the  eigenvalues 
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of  the  amplification  matrices  associated  with  each  of  the  three  parts 
must  be  less  than  or  equal  to  one  in  order  for  the  numerical  equations 
to  be  stable. 


Inviscid  Part 


The  inviscid  part  of  the  governing  equations  can  be  written  as 


|i+A|l+B3E  .0 

3t  3il 


The  amplification  matrix  of  this  equation,  as  given  by  MaeCormack 


(Ref  125),  is 


G - I 


- 14t  slna  + 


1 «.,*!  A-  a-c'15)  * L a-e't5)j  | ^ a-e1")  ♦ ^ 


where  cx  = k An,  £ = k AS,  and  1 is  the  unit  matrix.  The  von  Neumann 
l 2 


stability  condition  for  this  equation  is 


G < 1 


If  we  look  at  large  values  of  G,  then  the  G matrix  can  be  approximated 


G - I - iA'  - j (A’2  + B'2) 


where 


A1  * At  I * siaa  + 


k8in*] 


B'  - At 


[An  At,  J 
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If  a’  ia  the  eigenvalue  of  A'  and  b'  is  the  eigenvalue  of  B',  then  the 
stability  condition  is  satisfied  if 

1 (a’2  + b’2) 

4 b' 

Hence,  the  maximum  eigenvalue  of  the  matrix  A*  helps  determine  the 
maximum  allowable  time  increment  At.  Thus,  we  shall  let  sina  = sin0  = 1 

and  A'  « B ' . 

Richtm>er  and  Morton  (Ref  109)  introduced  a technique  by  which  the 
eigenvalues  of  A'  and  B'  may  be  found.  If  we  consider  an  axis  inclined 
at  an  angle  9 with  respect  to  the  n-axis  where  0 is  given  by 


1 


sin0 


_i 

An 


U) 


and  where  the  velocity  component  along  this  axis  is 


u'  *=  A cosO  + B sin0 

11  li 


then  the  A’  and  B'  matrices  can  be  written  as 


If  1 \2  . 

/ M 2 C 

/ ( An  ) 

I A?  J [ 

AcosO  + BsinO 


The  eigenvalues  for  this  matrix,  in  tensor  form,  arc 
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By  applying  the  MacCormack  scheme  to  the  linearized  form  of  this 
equation,  the  finite-difference  equation  becomes 
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Substitution  of  one  Fourier  component  of  the  solution 

U(n,C,  t)  = UQ(t)  exp  ^i(kj  n + k2  o] 

into  the  finite-difference  equation  gives  the  amplification  matrix, 
which  is  defined  as 
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where  a ■ k An.  8 ” k A£,  and  I is  the  unit  matrix.  If  we  let  c and  d 
l 2 

be  the  eigenvalues  of  the  C and  D matrices  respectively,  then  the 
von  Neumann  stability  criterion  becomes 


1 < 


1 

c + d 


The  maximum  eigenvalue  for  the  C matrix  is  the  larger  of  the  two  terms: 


where  8 is  the  maximum  local  value  of  the  normal  stress  damping 
function  in  the  x and  y directions.  For  the  I)  matrix,  the  maximum 
eigenvalue  is  the  larger  of  the  two  quantities: 


where  8 is  the  maximum  local  normal  stress  damping  in  the  x and  z 
2 

directions.  Thus,  the  maximum  allowable  time  increment  for  the  dif- 
fusion part  of  the  governing  equations  is  the  smaller  of  the  following 
two  terms: 


where  3 is  the  maximum  local  normal  stress  damping  vaiuu  in  the  three 

s 

Cartesian  coordinate  directions. 


Mixed  Derivative  Part 

The  mixed-derivative  part  of  the  governing  equations  is 
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By  applying  the  MacCormack  scheme  to  this  equation  results  in  the 
following  finite-difference  equation: 
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H *=  E + F 


and  where  the  coefficient  matrices  E and  F are  assumed  constant. 

If  we  substitute  the  Fourier  term 

u(n,£,t)  = Ufl(t)  exP  [l(kjn  + j 

into  the  finite-difference  equation,  the  amplification  matrix  becomes 
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whore  c is  an  adjustable  constant  loss  than  or  equal  to  one.  The 
magnitude  of  all  the  damping  term  constants  are  set  so  as  to  make 
the  damping  effective  in  the  lnviscid  I low  region  but  not  to 
appreciably  change  the  boundary  layer  behavior  or  modify  the  Reynolds 
number.  This  llneari/.ed  stability  criterion  mav  not  insure  numerical 
stability  in  all  cases  (Ref  12b),  but  it  should  be  valid  for  the 
experimental  cases  considered. 
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APPENDIX  D 


i: 


Elements  of  the  Linearized 
Matrix  Equations 

This  appendix  contains  a list  of  matrix  elements  used  In  the 
linearized  governing  equations.  These  matrix  elements  were  used 
to  determine  the  stability  criteria  of  the  numerical  integration. 


Aii  “ A22  * A33  * Ai» 4 = A55 


A12 
Al3 
Ai  4 
A2  5 
A35 
A4  5 
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As  3 
A54 


£ 3n 

p cix 

£ 3£ 

p 3y 

£ 

p 3z 

U«V£>  pg- 


(l+2ye) 

<1+™  PsT 


In  + u9ri 

3y  ''r3z 


li  1 * B22  = B33  = B44  = B55 
Biz  - P3— 

n H 

Bl3  “ p57 


Bl4 


3f. 

'3z 


n . E 3f. 
B2S  P ^ 


3f.  A 3f.  , 3£ 

J3x  + v3T  + W3z" 


200 




r 


I 


203 


APPENDIX  F. 


Subroutines  of  the  J) El ,TA  Code 

This  appendix  provides  a description  of  the  important  operations 
performed  by  each  subroutine  of  the  numerical  code  DELTA.  These  sub- 
routines are  also  used  in  the  DELTA1  and  CONE  programs.  The  order  ot 
the  subroutine  descriptions  follows  the  same  order  used  to  call  them 
in  the  DELTA  program. 

COORD 

Subroutine  COORD  calculates  the  first-order  derivatives  of  the 
coordinate  transformations  (as  defined  in  Appendix  A).  The  analytical 
values  of  these  derivatives  are  computed  for  every  grid  point  in  the 
computational  domain.  These  values  are  stored  in  common  block  arrays 
for  easy  access  during  the  numerical  integration. 

PREDICT 

The  next  subroutine  used  in  the  DELTA  program  is  PREDICT.  This 

subroutine  calculates  the  predictor  term  of  the  MacCormaek  finite- 

difference  scheme.  It  uses  the  ROUND  subroutine  to  compute  the 

boundary  conditions  on  the  surface  of  the  body  and  at  the  leading  edge 

of  the  delta  wing.  A double  DO  loop  is  entered  wherein  the  U " ' 

' » .1 

vector  is  calculated  for  all  computational  grid  points  except  the  free 
stream  boundary  points,  the  surface  grid  points,  and  the  grid  points 
opposite  the  plane  of  symmetry.  Within  this  DO  loop,  the  subroutines 
DAMPF,  PAM  PC. , VECTOR,  DECODE,  and  SOLVE  are  called  in  sequential  order 
to  perform  the  numerical  integration.  As  a result  of  this  partial 
integration,  the  flow  quantities  are  computed  tor  an  intermediate  time 
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BOUND 


Subroutine  BOUND  calculates  the  pressure  and  density  on  the  surtac 
of  the  bodv  and  at  the  leading  edge.  The  velocity  at  these  surface  gri 
points  Is  zero  and  the  surface  temperature  Is  provided  as  part  ot  the 
Input  data.  The  normal  momentum  equation  is  used  to  calculate  t lie 
surface  pressure  and  then  the  equation  of  state  is  used  to  determine 
the  density.  In  the  DELTA  program,  the  pressure  and  density  are  com- 
puted for  both  the  upper  and  lower  wing  surfaces,  as  described  in 
Chapter  3.  Second-order  accurate,  one-sided  forward  differences  are 
used  to  model  the  flow  gradients  normal  to  the  body  surface,  while 
second-order  central  differences  are  used  for  all  other  gradients. 

At  th  • leading  edge,  the  pressure  and  density  are  triple  valued. 
The  pressure  is  determined  for  the  upper,  side,  and  lower  surtaces  ot 
the  wing  bv  using  the  appropriate  normal  momentum  equal  ton.  I'ho  three 
values  of  density  are  computed  from  the  equation  of  state.  When  the 
numerical  Integration  is  applied  below  the  wing,  the  leading  edge  is 
represented  bv  lower  surface  values  ot  pressure  and  density.  Similarly 
when  the  numerical  Integration  occurs  above  the  wing,  the  upper  surface 
values  model  the  leading  edge  on  the  p-axis.  In  the  Dl'lTA l program, 
the  leading  edge  is  represented  bv  upper  surface  values  only. 

DANTE  and  PAMl'C 

DAM TV  and  DAMDO  subroutines  calculate  the  density  damping  terms 
for  the  n and  directions  respectively.  Vhese  terms  .lie  addl'd  togethe 
in  the  DECODE  subroutine  to  form  the  tourth-order  damping  terms  D^n 
(predictor  step  only)  and  d"  j (corrector  step  only!.  The  damp  lug 

l t .1 


coefficients,^  and  Cj  , are  read  into  the  program  as  part  of  the  input 
data.  The  damping  terms  normal  to  the  free  stream  boundary  surfaces 
and  symmetry  plane  and  within  two  grid  points  of  these  surfaces  are 
zero.  The  damping  at  the  surface  grid  points  is  also  zero.  All  the 
damping  term  values  are  stored  in  common  block  arrays  for  use  in  the 
DECODE  subroutine. 


VECTOR 

The  next  subroutine  used  in  the  DELTA  program  is  VECTOR.  This  sub- 
routine computes  the  values  of  the  F,G  and  H matrices  as  well  as  the 
heat  flux  and  shear  stress  terms  in  these  matrices.  The  finite-differ- 
ence quotients  used  to  model  and  q are  described  in  Chapter  3.  The 
second-order  derivatives  of  the  coordinate  transformations  (used  in  the 
H matrix)  are  calculated  by  calling  subroutine  COORDX.  The  values  of 
these  matrices  are  stored  in  common  block  arrays  for  use  in  the  DECODE 
subroutine . 

Subroutine  VECTOR  is  used  in  both  the  predictor  and  corrector 
steps.  By  changing  the  input  parameters  when  this  subroutine  is  called, 
the  finite-difference  quotients  can  be  easily  converted  for  use  in 
either  the  predictor  or  corrector  step.  These  input  parameters  are  also 
used  in  the  DECODE  subroutine. 


SOLVE 

Subroutine  SOLVE  is  used  to  solve  for  the  flow  quantities  in  the 

n+T  n+1 

U (predictor  step  only)  and  U . (corrector  step  only)  vectors, 

i * J i > J 

The  primitive  variables  u,v,w,p,  and  e are  determined  by  solving  the 


relationships  in  Eq  70. 


The  flow  quantities  are  then  stored  in  two- 


dimensional  arrays  for  use  in  the  next  integration  step. 
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SYM 
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Subroutine  SYM  uses  a DO  loop  to  calculate  the  primitive  flow  vari- 
ables on  the  opposite  side  of  the  plane  of  symmetry.  By  applying  Eqs 
79  and  80,  the  surface  boundary  conditions  and  the  flow  field  conditions 
are  determined  for  the  mirror  image  plane.  These  flow  quantities  are 
then  stored  in  common  block  arrays  for  use  in  the  next  integration  step. 

CORRECT 

The  next  subroutine  used  in  the  DELTA  program  is  CORRECT.  This 
subroutine  calculates  the  corrector  term  of  the  MacCormack  finite- 
difference  scheme.  It  uses  the  same  subroutines  as  are  used  in  the 
PREDICT  subroutine  in  the  same  sequential  order.  By  applying  different 
input  parameters,  this  subroutine  is  able  to  calculate  the  final  flow 
quantities  at  the  new  time  step  tn+^.  These  newly  calculated  t low 
parameters  are  then  stored  in  the  same  two-dimensional  arrays  ns  the 
old  quantities. 


APPENDIX  F 


Conical  Velocity  Compoin >nts 
and  Streamline  Plots 


This  appendix  presents  a detailed  discussion  on  the  conical  cross- 
flow  velocity  vectors  and  on  the  method,  used  by  Ghia  (Ref  139)  to 
determine  the  conical  streamline  contours  in  the  £-r|  plane.  The  com- 
puter software  for  the  streamline  plots  was  developed  by  Cooper  (Ref 
140)  . 


Conical  Velocity  Vectors 

Experimental  evidence  has  shown  that  the  supersonic  and  hypersonic 
flow  around  thin  delta  wings  is  nearly  conical  in  the  weak  interaction 
region.  In  delineating  the  conical  aspects  of  this  flow  field,  the 
velocity  components  in  the  spherical  coordinate  system  are  used  for 
flow  analysis.  ■ The  magnitude  of  the  spherical  cross-flow  velocity 
vectors  is 


+ ^ <F 
while  the  angle  of  the  flow  in  the  £~n  plane,  as  determined  by  Cooper 
and  Hankey  (Ref  142) , is 


tan  y 
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tan  (u)-if>) 


d£ 

dn 


0'2) 


where 


tan  a) 


UJ1 

cos  0 


When  these  values  are  transformed  into  the  conical  coordinate  system 
(delta  wing),  the  magnitude  and  direction  of  the  velocity  vectors 
become 
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f . .(“.n-v) 7 1 ** 

c L l + nJ  ■*■  CJ  J 


tan  y 


- I \ 

\ un-v  / 


The  steady  state  velocity  components  In  the  r)  and  C directions  are 


V - V cos  y 
H c 'c 


Vr  - V sin  Y 
£ c 'c 


Conical  Streamline  Plot 


The  cross-flow  velocity  vectors  projected  on  the  £-r)  plane  trace 
out  "pseudo"  conical  streamlines.  These  streamlines  represent  steadv 
state  path  lines  of  Incoming  fluid  particles.  These  path  lines  are 
determined  by  solving  the  differential  equations 


-I3  - V (t],t) 
dt  n 


mV  (n  n 
dt  V£  J 


The  modi  lied  Ku ler-Cauchy  method  (Kef  Isl)  is  used  to  solve  Equations 
E8  and  V'9  In  the  form 
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where  n and  £ are  the  projected  locations  of  the  next  point  on  the 
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and  where 
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Bv  averagim?  the  appropriate  velocity  components  at  (1  ,j  ) with 

* ★ 

those  at  (i,j)»  an  improved  estimate  of  n and  £ is  obtained.  This 
improved  estimate  is 


* 

n - n 


At 


At 


(F17) 


(FI  8) 


A A 

The  velocities  at  the  new  (l  ,J  ) are  recomputed,  by  linear  inter- 
polation, and  then  the  new  location  and  velocity  vector  are  used  as 
initial  conditions  for  calculating  tho  next  point.  This  integration 
process  is  repeated  until  the  maximum  number  of  time  steps  is  reached. 
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